1.Objetivos

  1. Obtendremos los datos reales de pacientes con cancer de mama de The Cancer Genome Atlas y prepararemos estos(convertirlos e imputarlos) para su consecuente análisis. Una vez los datos tengan el formato correcto, los estratificaremos de la forma en la que suelen tratarse en el ámbito médico. Por grupos tumorales, Luminal, grupos nodales, estado nodal y grupos de edad. 2

  2. Haremos un barrido rápido por algunas las opciones de análisis (estadístico 3.1, paramétrico 3.3 y no paramétrico 3.2) que podemos aplicar en análisis de supervivencia. Para ello usaremos principalmente el paquete survminer.

  3. Cuando tengamos el dataset final adecuado para trabajar con los datos correctos y las estratificaciones incluidas pasaremos a ver que variables nos interesan y cuales no, es decir, cuales son útiles para explicar los riesgos y la supervivencia. Este proceso lo llevaremos a cabo gracias al análisis multivariante del modelo semiparamétrico Cox. Realizaremos una búsqueda exhaustiva de las mejores covariantes (de forma univariante y multivariante), y después de descartar aquellos modelos que no cumplan la proporcionalidad de riesgos, nos quedaremos con aquella combinación de variables que explique mejor el riesgo. 3.4

  4. Una vez tengamos las variables con las que trabajaremos más en detalle, haremos un análisis exploratorio visual de como interactuan y explican el riesgo diferentes combinaciones de variables 3.5. De aquí sacaremos conclusiones de las cuales le haremos la consiguiente prueba de hipótesis 3.6.

  5. Analizaremos más en detalles los supuestos aceptados con más significancia. 3.7

2.Preparación datos

Instalación de paquetes y carga de librerias:

install.packages("survival")
install.packages("KMsurv")
install.packages("survMisc")
install.packages("survminer")
install.packages("UCSCXenaTools")
install.packages("TCGAretriever")
library(KMsurv)
library(survMisc)
library(survminer)
library(flexsurv)
library(dplyr)
library(UCSCXenaTools)
library(TCGAretriever) 
library(ggfortify)

Los extraeremos de TCGA, y después de imputarlos, los pondremos en el formato más conveniente para trabajarlos.

2.1. Importar los datos de TCGA

Obtendremos los datos de TCGA The Cancer Genome Atlas, usando el paquete TCGAretriever.

all_studies<-get_cancer_studies() #identifier, tittle and description of the study
str(all_studies)
## 'data.frame':    233 obs. of  3 variables:
##  $ cancer_study_id: chr  "paac_jhu_2014" "all_stjude_2016" "laml_tcga_pub" "laml_tcga_pan_can_atlas_2018" ...
##  $ name           : chr  "Acinar Cell Carcinoma of the Pancreas (Johns Hopkins, J Pathol 2014)" "Acute Lymphoblastic Leukemia (St Jude, Nat Genet 2016)" "Acute Myeloid Leukemia (TCGA, NEJM 2013)" "Acute Myeloid Leukemia (TCGA, PanCancer Atlas)" ...
##  $ description    : chr  "Mutation data from whole exome sequencing of 23 surgically resected pancreatic carcinomas with acinar differentiation." "Whole-genome and/or whole-exome sequencing was performed for ERG-altered B-ALL cases." "TCGA Acute Myeloid Leukemia, analysis of 200 adult cases; raw data at the <A HREF=\"https://tcga-data.nci.nih.g"| __truncated__ "Comprehensive TCGA PanCanAtlas (https://gdc.cancer.gov/about-data/publications/pancanatlas) data from 11k cases"| __truncated__ ...
head(all_studies[,2],20) #tittles of first 20 studies
##  [1] "Acinar Cell Carcinoma of the Pancreas (Johns Hopkins, J Pathol 2014)"  
##  [2] "Acute Lymphoblastic Leukemia (St Jude, Nat Genet 2016)"                
##  [3] "Acute Myeloid Leukemia (TCGA, NEJM 2013)"                              
##  [4] "Acute Myeloid Leukemia (TCGA, PanCancer Atlas)"                        
##  [5] "Acute Myeloid Leukemia (TCGA, Provisional)"                            
##  [6] "Adenoid Cystic Carcinoma (FMI, Am J Surg Pathl. 2014)"                 
##  [7] "Adenoid Cystic Carcinoma (MDA, Clin Cancer Res 2015)"                  
##  [8] "Adenoid Cystic Carcinoma (MSKCC, Nat Genet 2013)"                      
##  [9] "Adenoid Cystic Carcinoma (Sanger/MDA, JCI 2013)"                       
## [10] "Adenoid Cystic Carcinoma of the Breast (MSKCC, J Pathol. 2015)"        
## [11] "Adrenocortical Carcinoma (TCGA, PanCancer Atlas)"                      
## [12] "Adrenocortical Carcinoma (TCGA, Provisional)"                          
## [13] "Ampullary Carcinoma (Baylor College of Medicine, Cell Reports 2016)"   
## [14] "Bladder Cancer (MSKCC, Eur Urol 2014)"                                 
## [15] "Bladder Cancer (MSKCC, JCO 2013)"                                      
## [16] "Bladder Cancer (TCGA, Cell 2017)"                                      
## [17] "Bladder Cancer, Plasmacytoid Variant (MSKCC, Nat Genet 2016)"          
## [18] "Bladder Urothelial Carcinoma (BGI, Nat Genet 2013)"                    
## [19] "Bladder Urothelial Carcinoma (Dana Farber & MSKCC, Cancer Discov 2014)"
## [20] "Bladder Urothelial Carcinoma (TCGA, Nature 2014)"
head(all_studies[,1],20) #id of first 20 studies
##  [1] "paac_jhu_2014"                "all_stjude_2016"             
##  [3] "laml_tcga_pub"                "laml_tcga_pan_can_atlas_2018"
##  [5] "laml_tcga"                    "acyc_fmi_2014"               
##  [7] "acyc_mda_2015"                "acyc_mskcc_2013"             
##  [9] "acyc_sanger_2013"             "acbc_mskcc_2015"             
## [11] "acc_tcga_pan_can_atlas_2018"  "acc_tcga"                    
## [13] "ampca_bcm_2016"               "blca_mskcc_solit_2014"       
## [15] "blca_mskcc_solit_2012"        "blca_tcga_pub_2017"          
## [17] "blca_plasmacytoid_mskcc_2016" "blca_bgi"                    
## [19] "blca_dfarber_mskcc_2014"      "blca_tcga_pub"
#filter rows from all_studies which name contains breast
df.breast<-dplyr::select(filter(all_studies, grepl('Breast', all_studies$name)),everything())
df.breast$cancer_study_id
##  [1] "acbc_mskcc_2015"              "brca_metabric"               
##  [3] "breast_msk_2018"              "bfn_duke_nus_2015"           
##  [5] "brca_bccrc"                   "brca_broad"                  
##  [7] "brca_sanger"                  "brca_tcga_pub2015"           
##  [9] "brca_tcga_pub"                "brca_tcga_pan_can_atlas_2018"
## [11] "brca_tcga"                    "brca_bccrc_xenograft_2014"   
## [13] "brca_mbcproject_wagle_2017"
#from all the studies, we choose brca_tcga->brca_tcga_all
all_case_lists<-get_case_lists("brca_tcga")
colnames(all_case_lists)
## [1] "case_list_id"          "case_list_name"        "case_list_description"
## [4] "cancer_study_id"       "case_ids"
all_case_lists[,1:2] #id and name
##                        case_list_id
## 1           brca_tcga_3way_complete
## 2               brca_tcga_sequenced
## 3         brca_tcga_methylation_all
## 4                     brca_tcga_all
## 5  brca_tcga_protein_quantification
## 6                     brca_tcga_cna
## 7        brca_tcga_methylation_hm27
## 8       brca_tcga_methylation_hm450
## 9                    brca_tcga_mrna
## 10        brca_tcga_rna_seq_v2_mrna
## 11                   brca_tcga_rppa
## 12                 brca_tcga_cnaseq
##                                       case_list_name
## 1                                All Complete Tumors
## 2                               All Sequenced Tumors
## 3            All tumor samples with methylation data
## 4                                         All Tumors
## 5                 Protein Quantification (Mass Spec)
## 6                        Tumor Samples with CNA data
## 7         Tumor Samples with methylation data (HM27)
## 8        Tumor Samples with methylation data (HM450)
## 9  Tumor Samples with mRNA data (Agilent microarray)
## 10         Tumor Samples with mRNA data (RNA Seq V2)
## 11                      Tumor Samples with RPPA data
## 12        Tumor Samples with sequencing and CNA data
df.tcga<-get_clinical_data("brca_tcga_all")
str(df.tcga)
## 'data.frame':    1105 obs. of  110 variables:
##  $ CASE_ID                                   : chr  "TCGA-A7-A3J0-01" "TCGA-OL-A66N-01" "TCGA-AQ-A0Y5-01" "TCGA-E9-A22H-01" ...
##  $ AGE                                       : chr  "62" "59" "70" "42" ...
##  $ AJCC_METASTASIS_PATHOLOGIC_PM             : chr  "M0" "MX" "MX" "M0" ...
##  $ AJCC_NODES_PATHOLOGIC_PN                  : chr  "N0" "N3" "N2a" "N1" ...
##  $ AJCC_PATHOLOGIC_TUMOR_STAGE               : chr  "Stage IIA" "Stage IIIC" "Stage IIIA" "Stage IIB" ...
##  $ AJCC_STAGING_EDITION                      : chr  "7th" "7th" "7th" "7th" ...
##  $ AJCC_TUMOR_PATHOLOGIC_PT                  : chr  "T2" "T3" "T2" "T2" ...
##  $ BRACHYTHERAPY_TOTAL_DOSE_POINT_A          : chr  "" "" "" "" ...
##  $ CANCER_TYPE                               : chr  "Breast Cancer" "Breast Cancer" "Breast Cancer" "Breast Cancer" ...
##  $ CANCER_TYPE_DETAILED                      : chr  "Breast Invasive Mixed Mucinous Carcinoma" "Breast Invasive Lobular Carcinoma" "Breast Invasive Ductal Carcinoma" "Breast Invasive Ductal Carcinoma" ...
##  $ CENT17_COPY_NUMBER                        : chr  "" "" "" "" ...
##  $ DAYS_TO_BIRTH                             : chr  "-22672" "-21608" "-25793" "-15504" ...
##  $ DAYS_TO_COLLECTION                        : chr  "78" "505" "212" "38" ...
##  $ DAYS_TO_DEATH                             : chr  "" "" "172" "" ...
##  $ DAYS_TO_INITIAL_PATHOLOGIC_DIAGNOSIS      : chr  "0" "0" "0" "0" ...
##  $ DAYS_TO_LAST_FOLLOWUP                     : chr  "76" "413" "" "0" ...
##  $ DFS_MONTHS                                : chr  "10.28" "26.02" "" "40.47" ...
##  $ DFS_STATUS                                : chr  "DiseaseFree" "DiseaseFree" "" "DiseaseFree" ...
##  $ DISEASE_CODE                              : chr  "" "" "" "" ...
##  $ ER_POSITIVITY_SCALE_OTHER                 : chr  "" "" "" "" ...
##  $ ER_POSITIVITY_SCALE_USED                  : chr  "" "" "" "" ...
##  $ ER_STATUS_BY_IHC                          : chr  "Positive" "Positive" "Positive" "Positive" ...
##  $ ER_STATUS_IHC_PERCENT_POSITIVE            : chr  "90-99%" "80-89%" "70-79%" "<10%" ...
##  $ ETHNICITY                                 : chr  "NOT HISPANIC OR LATINO" "NOT HISPANIC OR LATINO" "NOT HISPANIC OR LATINO" "NOT HISPANIC OR LATINO" ...
##  $ FIRST_SURGICAL_PROCEDURE_OTHER            : chr  "" "" "" "" ...
##  $ FORM_COMPLETION_DATE                      : chr  "3/13/12" "6/11/13" "4/7/11" "8/2/11" ...
##  $ FRACTION_GENOME_ALTERED                   : chr  "0.224023880" "0.209382668" "0.496072001" "0.103942673" ...
##  $ HER2_AND_CENT17_CELLS_COUNT               : chr  "" "" "" "" ...
##  $ HER2_AND_CENT17_SCALE_OTHER               : chr  "" "" "" "" ...
##  $ HER2_CENT17_RATIO                         : chr  "" "1.1" "3.8" "" ...
##  $ HER2_COPY_NUMBER                          : chr  "" "" "" "" ...
##  $ HER2_FISH_METHOD                          : chr  "" "" "" "" ...
##  $ HER2_FISH_STATUS                          : chr  "" "Negative" "Positive" "" ...
##  $ HER2_IHC_PERCENT_POSITIVE                 : chr  "" "" "" "<10%" ...
##  $ HER2_IHC_SCORE                            : chr  "1" "" "2" "" ...
##  $ HER2_POSITIVITY_METHOD_TEXT               : chr  "" "" "" "" ...
##  $ HER2_POSITIVITY_SCALE_OTHER               : chr  "" "" "" "" ...
##  $ HISTOLOGICAL_DIAGNOSIS                    : chr  "Mucinous Carcinoma" "Infiltrating Lobular Carcinoma" "Infiltrating Ductal Carcinoma" "Infiltrating Ductal Carcinoma" ...
##  $ HISTOLOGICAL_SUBTYPE                      : chr  "" "" "" "" ...
##  $ HISTORY_NEOADJUVANT_TRTYN                 : chr  "No" "No" "Yes" "No" ...
##  $ HISTORY_OTHER_MALIGNANCY                  : chr  "No" "No" "Yes" "No" ...
##  $ ICD_10                                    : chr  "C50.9" "C50.9" "C50.9" "C50.9" ...
##  $ ICD_O_3_HISTOLOGY                         : chr  "8480/3" "8520/3" "8500/3" "8500/3" ...
##  $ ICD_O_3_SITE                              : chr  "C50.9" "C50.9" "C50.9" "C50.9" ...
##  $ IHC_HER2                                  : chr  "Negative" "" "Positive" "Positive" ...
##  $ IHC_SCORE                                 : chr  "" "" "" "" ...
##  $ INFORMED_CONSENT_VERIFIED                 : chr  "YES" "YES" "YES" "YES" ...
##  $ INITIAL_PATHOLOGIC_DX_YEAR                : chr  "2011" "2011" "2010" "2011" ...
##  $ INITIAL_WEIGHT                            : chr  "110" "310" "160" "500" ...
##  $ IS_FFPE                                   : chr  "NO" "NO" "NO" "NO" ...
##  $ LYMPH_NODES_EXAMINED                      : chr  "YES" "YES" "YES" "YES" ...
##  $ LYMPH_NODES_EXAMINED_HE_COUNT             : chr  "0" "13" "5" "1" ...
##  $ LYMPH_NODES_EXAMINED_IHC_COUNT            : chr  "" "" "" "0" ...
##  $ LYMPH_NODE_EXAMINED_COUNT                 : chr  "2" "15" "17" "8" ...
##  $ MARGIN_STATUS_REEXCISION                  : chr  "" "" "" "" ...
##  $ MENOPAUSE_STATUS                          : chr  "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" "Pre (<6 months since LMP AND no prior bilateral ovariectomy AND not on estrogen replacement)" ...
##  $ METASTATIC_SITE                           : chr  "" "" "" "" ...
##  $ METASTATIC_SITE_OTHER                     : chr  "" "" "" "" ...
##  $ METASTATIC_TUMOR_INDICATOR                : chr  "NO" "" "NO" "NO" ...
##  $ METHOD_OF_INITIAL_SAMPLE_PROCUREMENT      : chr  "Core needle biopsy" "Core needle biopsy" "Core needle biopsy" "Tumor resection" ...
##  $ METHOD_OF_INITIAL_SAMPLE_PROCUREMENT_OTHER: chr  "" "" "" "" ...
##  $ MICROMET_DETECTION_BY_IHC                 : chr  "NO" "" "NO" "YES" ...
##  $ MUTATION_COUNT                            : chr  "24" "" "35" "18" ...
##  $ NEW_TUMOR_EVENT_AFTER_INITIAL_TREATMENT   : chr  "" "NO" "" "" ...
##  $ NTE_CENT_17_HER2_RATIO                    : chr  "" "" "" "" ...
##  $ NTE_ER_IHC_INTENSITY_SCORE                : chr  "" "" "" "" ...
##  $ NTE_ER_STATUS                             : chr  "" "" "" "" ...
##  $ NTE_ER_STATUS_IHC_POSITIVE                : chr  "" "" "" "" ...
##  $ NTE_HER2_FISH_STATUS                      : chr  "" "" "" "" ...
##  $ NTE_HER2_POSITIVITY_IHC_SCORE             : chr  "" "" "" "" ...
##  $ NTE_HER2_STATUS                           : chr  "" "" "" "" ...
##  $ NTE_HER2_STATUS_IHC_POSITIVE              : chr  "" "" "" "" ...
##  $ NTE_PR_IHC_INTENSITY_SCORE                : chr  "" "" "" "" ...
##  $ NTE_PR_STATUS_BY_IHC                      : chr  "" "" "" "" ...
##  $ NTE_PR_STATUS_IHC_POSITIVE                : chr  "" "" "" "" ...
##  $ OCT_EMBEDDED                              : chr  "FALSE" "TRUE" "FALSE" "TRUE" ...
##  $ ONCOTREE_CODE                             : chr  "IMMC" "ILC" "IDC" "IDC" ...
##  $ OS_MONTHS                                 : chr  "10.28" "26.02" "5.65" "40.47" ...
##  $ OS_STATUS                                 : chr  "LIVING" "LIVING" "DECEASED" "LIVING" ...
##  $ OTHER_PATIENT_ID                          : chr  "DA70CF7E-0E61-4C72-B4C5-C408569D11B8" "E9002F4E-BA0B-4738-9368-DD9614A5DA8A" "d0b78f3f-a198-437a-ab8c-204345d3b75d" "161917b8-88ad-407b-ade6-d6b98478b359" ...
##  $ OTHER_SAMPLE_ID                           : chr  "025EAC6D-539E-4E43-A6F7-4A2D2DEF571D" "34E08E44-8C58-48BD-B682-1EDE7BAEAD3C" "f5675b71-ed1b-4fdb-be07-d627dc98ed88" "9ff0477c-10cd-47d6-8d0d-4e763ea185a8" ...
##  $ PATHOLOGY_REPORT_FILE_NAME                : chr  "TCGA-A7-A3J0.A5AED5F0-A144-4615-AF78-539B8523D15F.pdf" "TCGA-OL-A66N.CDB6A529-DE8B-48F6-A569-D07251D139C3.pdf" "TCGA-AQ-A0Y5.5BC01981-366E-44EC-A986-087DB7904E5C.pdf" "TCGA-E9-A22H.672E4D81-DA47-495C-A58B-BD2B4FEA2380.pdf" ...
##  $ PATHOLOGY_REPORT_UUID                     : chr  "A5AED5F0-A144-4615-AF78-539B8523D15F" "CDB6A529-DE8B-48F6-A569-D07251D139C3" "5BC01981-366E-44EC-A986-087DB7904E5C" "672E4D81-DA47-495C-A58B-BD2B4FEA2380" ...
##  $ PATH_MARGIN                               : chr  "Negative" "Negative" "Negative" "Negative" ...
##  $ PHARMACEUTICAL_TX_ADJUVANT                : chr  "" "NO" "" "" ...
##  $ PRIMARY_SITE                              : chr  "Right" "Left" "Right" "Left Upper Outer Quadrant" ...
##  $ PROJECT_CODE                              : chr  "" "" "" "" ...
##  $ PROSPECTIVE_COLLECTION                    : chr  "YES" "NO" "YES" "YES" ...
##  $ PR_POSITIVITY_DEFINE_METHOD               : chr  "" "" "" "" ...
##  $ PR_POSITIVITY_IHC_INTENSITY_SCORE         : chr  "" "" "" "" ...
##  $ PR_POSITIVITY_SCALE_OTHER                 : chr  "" "" "" "" ...
##  $ PR_POSITIVITY_SCALE_USED                  : chr  "" "" "" "" ...
##  $ PR_STATUS_BY_IHC                          : chr  "Positive" "Negative" "Positive" "Positive" ...
##  $ PR_STATUS_IHC_PERCENT_POSITIVE            : chr  "90-99%" "" "40-49%" "10-19%" ...
##  $ RACE                                      : chr  "WHITE" "WHITE" "WHITE" "WHITE" ...
##  $ RADIATION_TREATMENT_ADJUVANT              : chr  "" "YES" "" "" ...
##  $ RETROSPECTIVE_COLLECTION                  : chr  "NO" "YES" "NO" "NO" ...
##  $ SAMPLE_TYPE                               : chr  "Primary" "Primary" "Primary" "Primary" ...
##  $ SAMPLE_TYPE_ID                            : chr  "1" "1" "1" "1" ...
##   [list output truncated]

Después de una vista preliminar de todos los datos que contiene nos quedamos con los que más nos interesan (explicados con detalle más adelante) y los guardamos en un dataframe.

str(df.tcga$AGE) #age
##  chr [1:1105] "62" "59" "70" "42" "69" "61" "50" "80" "59" "65" "79" ...
str(df.tcga$MENOPAUSE_STATUS) #menostat
##  chr [1:1105] "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" ...
str(df.tcga$AJCC_PATHOLOGIC_TUMOR_STAGE) #tgrade
##  chr [1:1105] "Stage IIA" "Stage IIIC" "Stage IIIA" "Stage IIB" ...
str(df.tcga$LYMPH_NODE_EXAMINED_COUNT)  #pnodes
##  chr [1:1105] "2" "15" "17" "8" "4" "4" "30" "10" "" "4" "2" "10" "12" ...
str(df.tcga$ER_STATUS_BY_IHC) #strec
##  chr [1:1105] "Positive" "Positive" "Positive" "Positive" "Positive" ...
str(df.tcga$PR_STATUS_BY_IHC) #progrec
##  chr [1:1105] "Positive" "Negative" "Positive" "Positive" "Positive" ...
str(df.tcga$IHC_HER2) #her2
##  chr [1:1105] "Negative" "" "Positive" "Positive" "Equivocal" ...
str(df.tcga$OS_MONTHS) #time (months)
##  chr [1:1105] "10.28" "26.02" "5.65" "40.47" "24.47" "97.77" "85.58" ...
str(df.tcga$OS_STATUS) #cens
##  chr [1:1105] "LIVING" "LIVING" "DECEASED" "LIVING" "LIVING" "LIVING" ...
#Remove inteterminate and empty values 
  #remove os_months=0
df.tcga<-df.tcga[df.tcga$OS_MONTHS>0,]
  #remove menostat with values Indeterminate or empty
df.tcga<-df.tcga[as.factor(df.tcga$MENOPAUSE_STATUS)!="Indeterminate (neither Pre or Postmenopausal)",]
df.tcga<-df.tcga[as.factor(df.tcga$MENOPAUSE_STATUS)!="",]
  #remove tgrade empty
df.tcga<-df.tcga[as.factor(df.tcga$AJCC_PATHOLOGIC_TUMOR_STAGE)!="",]
  #remove HER, PR, ER empty or indeterminate
df.tcga<-df.tcga[as.factor(df.tcga$ER_STATUS_BY_IHC)!="",]
df.tcga<-df.tcga[as.factor(df.tcga$PR_STATUS_BY_IHC)!="",]
df.tcga<-df.tcga[as.factor(df.tcga$PR_STATUS_BY_IHC)!="Indeterminate",]
df.tcga<-df.tcga[as.factor(df.tcga$IHC_HER2)!="",]
df.tcga<-df.tcga[as.factor(df.tcga$IHC_HER2)!="Indeterminate",]

#Create dataframe
df.tcga$OS_MONTHS<-as.numeric(df.tcga$OS_MONTHS)
df.final<-data.frame(df.tcga$AGE,df.tcga$MENOPAUSE_STATUS,df.tcga$AJCC_PATHOLOGIC_TUMOR_STAGE,df.tcga$LYMPH_NODES_EXAMINED_HE_COUNT,df.tcga$ER_STATUS_BY_IHC,df.tcga$PR_STATUS_BY_IHC,df.tcga$IHC_HER2,df.tcga$OS_MONTHS,df.tcga$OS_STATUS)
#name dataframe columns
names(df.final)<-c("age", "menostat",  "tgrade", "pnodes", "ER", "PR", "HER2", "time", "cens")

2.2. Prepararlos para su análisis (Conversión e imputación)

Dado que muchos de los datos no están en un formato válido “”, NA,…, los imputaremos.

Para los datos numéricos, en este caso edad y número de nodos sustituiremos los valores perdidos con la media del resto.

Para los datos categóricos haremos dos secciones, aquellos que serán imputados con la moda (menostat y tgrade), y aquellos que en caso de ausencia serán interpretados como negativos (ER, HER2, PR). El valor de censura será 0 en caso de no aparecer.

impute.cols<-function(df) {
  
  #NUMERIC->MEAN
  imputar.numericos<-function(df, col) {
    mean.value<-mean(as.integer(na.omit(df[,col])))
    df[,col][df[,col]==""]<-floor(mean.value)
    df[,col]<-as.numeric(as.character(df[,col])) 
    return(df[,col])
  }

  df$age<-imputar.numericos(df,"age")
  df$pnodes<-imputar.numericos(df,"pnodes")

  #CATEGORICALS
  cat("LEVELS OF CATEGORICALS: \n")
  #menostat
  mode.menostat<-names(table(df$menostat)[which.max(table(df$menostat))])
  mode.menostat
  cat("Menostat levels: ", levels(df$menostat), "\n") # Peri as (Post)
  df$menostat<-factor(df$menostat, labels=c("Post", "Post", "Pre"))
  #tgrade
  mode.tgrade<-names(table(df$tgrade)[which.max(table(df$tgrade))])
  mode.tgrade
  cat("Tumor grade levels: ", levels(df$tgrade), "\n") #X as mode (II)
  df$tgrade<-factor(df$tgrade, labels=c("I","I", "I", "II", "II", "II", "III", "III", "III", "III", "IV", "II"))
  #cens 
  cat("cens levels: ", levels(df.final$cens), "\n")
  df$cens<-as.numeric(factor(df$cens, labels=c("1", "0"))) -1 #factor starts at 0
  #ER
  cat("Estrogen Receptor levels: ", levels(df$ER), "\n") 
  df$ER<-factor(df$ER, labels=c("Negative", "Positive"))
  #HER2
  cat("HER2 levels: ",levels(df$HER2), "\n") #Equivocal as Negative
  df$HER2<-factor(df$HER2, labels=c("Negative", "Negative", "Positive"))
  #PR
  cat("Progesterone levels: ", levels(df$PR), "\n")

df$PR<-factor(df$PR, labels=c("Negative", "Positive"))
  
  return(df)
}

df.final.imp<-impute.cols(df.final)
## LEVELS OF CATEGORICALS: 
## Menostat levels:  Peri (6-12 months since last menstrual period) Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy) Pre (<6 months since LMP AND no prior bilateral ovariectomy AND not on estrogen replacement) 
## Tumor grade levels:  Stage I Stage IA Stage IB Stage II Stage IIA Stage IIB Stage III Stage IIIA Stage IIIB Stage IIIC Stage IV Stage X 
## cens levels:  DECEASED LIVING 
## Estrogen Receptor levels:  Negative Positive 
## HER2 levels:  Equivocal Negative Positive 
## Progesterone levels:  Negative Positive
df.final.imp[1:20,]
##    age menostat tgrade pnodes       ER       PR     HER2   time cens
## 1   62     Post     II      0 Positive Positive Negative  10.28    1
## 2   70     Post    III      5 Positive Positive Positive   5.65    0
## 3   42      Pre     II      1 Positive Positive Positive  40.47    1
## 4   69     Post      I      0 Positive Positive Negative  24.47    1
## 5   61     Post      I      0 Positive Positive Negative  97.77    1
## 6   50      Pre    III      9 Positive Positive Positive  85.58    1
## 7   80     Post     II      2 Positive Positive Positive  12.19    1
## 8   59     Post    III      6 Positive Positive Negative  32.98    1
## 9   65     Post     II      0 Positive Positive Negative 160.78    1
## 10  49     Post     II      2 Positive Negative Negative  14.78    1
## 11  62     Post     II      6 Positive Negative Negative  10.58    1
## 12  64     Post     II      0 Positive Positive Negative  50.82    1
## 13  64     Post     II      3 Positive Positive Negative  56.90    1
## 14  58     Post     II      0 Negative Positive Positive  36.17    1
## 15  85     Post     II      6 Positive Positive Negative  45.04    1
## 16  79     Post     II      1 Positive Positive Negative  44.78    1
## 17  54     Post    III      9 Positive Positive Negative  67.81    1
## 18  36      Pre    III      4 Positive Positive Positive  17.05    1
## 19  53     Post     II      2 Negative Negative Negative  39.52    1
## 20  49      Pre     II      0 Positive Positive Negative  13.90    1
BC<-df.final.imp

2.3.Datos finales

  • age: Edad de la paciente en años.
#Living patients
qplot(BC$age, geom="histogram")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Deceased patients
BC.age.nocens<-BC[BC$cens==1,]$age
qplot(BC.age.nocens, geom="histogram")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  • menostat: Estado de la menopausia, premenopausal o postmenopausal.
table(BC$menostat)
## 
## Post  Pre 
##  621  196

Se encuentran en mayor proporción las pacientes Postmenopausicas.

  • tgrade: Grado del tumor I, II, III, IV.

Se determina por el valor de TNM (Tumor tamaño y localización. Node ¿Se ha expandido a nodos linfáticos? ¿Donde y cuantos?. Metastasis ¿Ha habido metástasis a otras partes del cuerpo?).

Siendo:

  • Grado I: 20 millimeters (mm) o menor en su area más ancha.

  • Grado II: El tamaño del tumor es mayor a 20mm y menor a 50mm.

  • Grado III: Mayor a 50 mm

  • Grado IV: Puede ser alguno de estos casos: ha crecido en la pared del pecho(TIVa), en la piel(TIVb) o ambos(TIVc), o el cancer es inflamatorio (TIVd).

table(BC$tgrade)
## 
##   I  II III  IV 
## 142 484 179  12
  • pnodes: Número de nodos linfáticos encontrados.

  • progrec: Receptor de progesterona (hormona esteroide C-21 involucrada en el ciclo menstrual, el embarazoy la embriogénesis) en \([fmol]\) (\(10e-15\) moles)

  • estrec: Receptor de estrógeno (hormonas sexuales esteroideas).

Estrógeno y Progesterona en el cancer de mama: Si las células del tumor tienen receptores de estrógeno, se dice que el cáncer es positivo en cuanto a receptores de estrógeno, sensible al estrógeno o que responde al estrógeno (ER-positivo), si no tienen, son ER-negativos y no usan el estrógeno para crecer. De forma semejante, si las células del tumor tienen receptores de progesterona, se dice que el cáncer es positivo en cuanto a receptores de progesterona (PR-positivo o PgR-positivo). Aproximadamente 80 % de los cánceres de mama tienen receptores de estrógeno (La mayoría de los cánceres de seno ER-positivos son también PR-positivos)

Si no tienen ninguno de los dos son HR negativos, receptores negativos de hormonas.

table(BC$Luminal)
## < table of extent 0 >

Aunque no es mucha la diferencia, podemos ver que efectivamente si hay más casos que presentan estrógeno de los que presentan progesterona.

  • time: Tiempo transcurrido en estudio, en meses.

Si visualizamos el tiempo en un histograma, podremos ver claramente que no tiene una distribución (a diferencia de por ejemplo, la edad como hemos visto anteriormente), por lo que tendremos que aplicar métodos no paramétricos.

qplot(BC$time, geom="histogram")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  • cens: Si es un caso censurado(0) o no(1), es decir, si no sabemos el momento exacto de ocurrencia del evento, ya sea por abandono del estudio o porque no este no se produzca.

2.3.1 Estratificaciones de variables

2.3.1.1.Categóricas

Luminal:

Ya hemos explicado al incio de este apartado los significados de los grupos de hormonas ER, HER2 y PR.

  • LuminalA: ER positivo y HER2 negativo. El grado tumoral es 1 o 2. Entre el 30-70% de los tumores de cancer de mama son de este tipo, y tienen a tener el mejor pronóstico.

-Luminal B: ER positivo. Lo tienen un 10-20% de las pacientes.

-TN (Triple negative/basal-like): ER-negative y PR-negative. Lo tienen un 15-20%, suelen ser agresivos y tener peor diagnóstico que los anteriores.

#LUMINAL
make.luminal<-function(df) {
  vectorluminals <-list()
  for(i in 1:nrow(df)) {
    hormone_receptor_positive=(df[i,]$ER=="Positive" || df[i,]$PR=="Positive")
    hormone_receptor_negative=(df[i,]$ER=="Negative" && df[i,]$PR=="Negative")
    if(hormone_receptor_positive && df[i,]$HER2=="Negative") vectorluminals[i]<-"LuminalA"
    else if(hormone_receptor_positive)  vectorluminals[i]<-"LuminalB" 
    else if(df[i,]$ER=="Negative" && df[i,]$PR=="Negative" && df[i,]$HER2=="Negative") vectorluminals[i]<-"TN"
    else if(hormone_receptor_negative && df[i,]$HER2=="Positive") vectorluminals[i]<-"HER2-enriched"
    else if(hormone_receptor_positive && df[i,]$HER2=="Negative") vectorluminals[i]<-"Normal-line"
  }
  df$Luminal<-as.factor(unlist(vectorluminals))
  return(df)
}

BC<-make.luminal(BC)


#NODAL AND NODE STATUS
make.nodalStatus<-function(df) {
#five groups: 0, negative; 1, one positive; 2, two to three positive; 3, four to nine positive; and 4, >9 nine positive
  vectorGroups<-list()
  nodeSt<-list()
  for(n in 1:length(df$pnodes)) {
    numnod<-df$pnodes[n]
    if(numnod==0) {
      vectorGroups[n]="G0"
      nodeSt[n]="Node negative"
    } 
    else if(numnod>0) {
      nodeSt[n]="Node positive"
      if(numnod==1) vectorGroups[n]="G1"
      else if(numnod==2 || numnod==3) vectorGroups[n]="G2"
      else if(numnod>=4 &&numnod<=9) vectorGroups[n]="G3"
      else if(numnod>9) vectorGroups[n]="G4"
    }
  }
  df$NodalStatus<-as.factor(unlist(vectorGroups))
  df$NodeStatus<-as.factor(unlist(nodeSt))
  return(df)
}

BC<-make.nodalStatus(BC)

#Otra agrupación distinta
max=max(BC$pnodes)
BC$NodalStatus2<-cut(BC$pnodes, breaks=c(0,1,3,5, max), include.lowest=T)

Veamos si la estratificación de Luminal cumple las proporciones esperadas:

prop.table(table(BC$Luminal))
## 
## HER2-enriched      LuminalA      LuminalB            TN 
##    0.04161567    0.65850673    0.13463892    0.16523868

Se adaptan perfectamente a los porcentajes descritos.

2.3.1.1.Numéricas

Estratificaremos la edad mediante clustering. También podriamos haber agrupado de esta forma el número de nodos, pero en su lugar nos hemos guiado por la definición medica con los grupos ya establecidos y lo hemos hecho de la misma manera.

Primero veremos como de bueno es el resultado para cada número de centroides, y procederemos a hacer los grupos de edad en una nueva variable.

library(factoextra)
## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
#Data
data<-BC
data$d1<-rep(1,length(BC$age))
#Convert factor variables to numeric 'dummy' variables
data.num=model.matrix(d1+age~age, data)
#número óptimo de clusters con su puntuación silohuette
fviz_nbclust(data.frame(data$age,data$d1), kmeans, method = "silhouette") +geom_vline(xintercept = 2, linetype = 2)

Vemos que el número óptimo de clusters para nuestra variable edad según el silhouette index es 2, 6 o 10, con un valor medio de 0.6. Aún así, haremos 3 clusters ya que queremos grupos más diferenciados y también dan un resultado de clusterización bueno.

#Assign clusters
data$cl= kmeans(data.num,3)$cluster

#df<-na.omit(data.frame(data$age,data$d1))
#fviz_cluster(kmeans(data.num,3), df, geom="point") + ggtitle("k = 3")

# Vemos como se han distribuido los datos
table(data$cl)
## 
##   1   2   3 
## 334 154 329
table(data$age[data$cl==1])
## 
## 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 
## 19 26 12 24 26 14 26 38 26 23 19 16 14 19 17 15
table(data$age[data$cl==2])
## 
## 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 
## 15 10 11 10  9 11  9 10 14 10  6  4  4  8  5  1  2  4  3  8
table(data$age[data$cl==3])
## 
## 26 27 28 29 30 31 32 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 
##  1  1  1  4  1  1  1  8  6  8  6  9  9 11 16 10  8  9 20 22 20 21 17 25 21 
## 52 53 54 
## 21 22 30
BC$GrupoEdadCluster<-cut(BC$age, breaks=c(25,50,66,90))
table(BC$GrupoEdadCluster)
## 
## (25,50] (50,66] (66,90] 
##     235     363     219
#Otra agrupación distinta
BC$GrupoEdad2<-cut(BC$age, breaks=c(25,37,50,90))

3. Análisis

Creamos un objeto supervivencia con Surv() a partir del tiempo hasta el evento y el indicador de censura (apareceran los tiempos y un \(+\) alado de los casos que sí estén censurados).

surv.obj<-Surv(BC$time, BC$cens) 
surv.obj[1:50]
##  [1]  10.28    5.65+  40.47   24.47   97.77   85.58   12.19   32.98 
##  [9] 160.78   14.78   10.58   50.82   56.90   36.17   45.04   44.78 
## [17]  67.81   17.05   39.52   13.90   31.34   54.17+   0.33   14.72 
## [25]  98.26   83.80+ 115.60   22.31   17.51   41.59   21.29    9.99 
## [33]  15.44   36.83   14.72   15.47   91.92+  16.03   30.98+   9.00 
## [41]  15.14   98.19   10.78   12.48   38.93   25.20   23.88   18.86 
## [49]  44.38   82.56

3.1. Análisis estadístico

3.1.1. Survfit

Con survfit calculamos el estimador que se le indique, como ya hemos visto que no es una distribución paramétrica en la variable tiempo, usaremos estimadores no paramétricos.

En survfit se indica objeto surv ~ nombre de las covariables o 1, en cuyo caso lo hará para la población global, usando todas las pacientes.

Kaplain-Meier, el que usaremos principalmente es un estimador no paramétrco que calcula la supervivencia (de los pacientes cuyo evento se ha observado), en proporciones exactas.

bc.fit.km<-survfit(surv.obj~1, data=BC, type="kaplan-meier")
bc.fit.km
## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "kaplan-meier")
## 
##       n  events  median 0.95LCL 0.95UCL 
##   817.0   738.0    27.8    24.5    31.0
bc.fit.fh<-survfit(surv.obj~1, data=BC, type="fleming-harrington")
bc.fit.fh
## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "fleming-harrington")
## 
##       n  events  median 0.95LCL 0.95UCL 
##   817.0   738.0    27.8    24.5    31.3

El el objeto de tipo survfit observamos el número de pacientes, la media de supervicia y el intervalo de confianza \(95%\).

Con summary() podemos ver los datos que contiene, lo haremos para algunos tiempos en concreto.

summary(bc.fit.km, time=c(50,200))
## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "kaplan-meier")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    50    222     543  0.30567 0.01665     0.274713       0.3401
##   200      1     195  0.00238 0.00235     0.000344       0.0165

time el tiempo de la observación, n.risk el número de sujetos en riesgo en ese tiempo, n.evento el número de sujetos que presentaron el evento, survival la estimación de la función de supervivencia (que es menor conforme aumenta el tiempo por lo que explicaremos más adelante). Y dado que estamos estimando la población global solo con una muestra, esto conlleva un arror asociado std.err, con sus correspondientes intervalos de cofidencia CI 95%.

También tenemos la opción de tenerlo como un dataframe usando la función fortity()

bc.fit.km.df<-fortify(bc.fit.km) 

Aprovechamos el dataframe para observar la clara diferencia de los primeros datos con los últimos (ordenados por tiempo ascendente). Como veremos y explicaremos más adelante en las gráficas, al principio la supervivencia es prácticamente 1, así como los intervalos, y un error muy bajo, valores que cambian inversamente llegando al final de los tiempos de estudio.

bc.fit.km.df[1:10,]
##    time n.risk n.event n.censor      surv     std.err     upper     lower
## 1  0.03    817       1        0 0.9987760 0.001224740 1.0000000 0.9963814
## 2  0.16    816       2        0 0.9963280 0.002123916 1.0000000 0.9921891
## 3  0.23    814       1        0 0.9951040 0.002453995 0.9999018 0.9903293
## 4  0.30    813       1        0 0.9938800 0.002745339 0.9992423 0.9885466
## 5  0.33    812      15        0 0.9755202 0.005542106 0.9861744 0.9649811
## 6  0.53    797       1        0 0.9742962 0.005682535 0.9852081 0.9635051
## 7  0.62    796       1        0 0.9730722 0.005819916 0.9842354 0.9620356
## 8  0.79    795       1        0 0.9718482 0.005954464 0.9832567 0.9605722
## 9  0.99    794       4        0 0.9669523 0.006467811 0.9792880 0.9547719
## 10 1.02    790       4        0 0.9620563 0.006947980 0.9752470 0.9490440
bc.fit.km.df[(nrow(bc.fit.km.df)-10):nrow(bc.fit.km.df),]
##       time n.risk n.event n.censor        surv   std.err      upper
## 630 125.07     11       1        0 0.023817389 0.2705758 0.04047703
## 631 129.99     10       1        0 0.021435650 0.2903831 0.03787139
## 632 131.57      9       1        0 0.019053911 0.3133867 0.03521595
## 633 132.95      8       1        0 0.016672172 0.3406881 0.03250771
## 634 133.11      7       1        0 0.014290434 0.3740026 0.02974383
## 635 134.03      6       1        0 0.011908695 0.4161866 0.02692296
## 636 134.30      5       1        0 0.009526956 0.4724524 0.02404952
## 637 136.63      4       1        0 0.007145217 0.5536647 0.02114928
## 638 140.44      3       1        0 0.004763478 0.6879035 0.01834293
## 639 160.78      2       1        0 0.002381739 0.9865147 0.01646709
## 640 216.59      1       0        1 0.002381739 0.9865147 0.01646709
##            lower
## 630 0.0140145656
## 631 0.0121328297
## 632 0.0103092926
## 633 0.0085506267
## 634 0.0068658427
## 635 0.0052675119
## 636 0.0037740005
## 637 0.0024139884
## 638 0.0012370282
## 639 0.0003444858
## 640 0.0003444858

También podemos ver como sería en el caso de estratificarse con alguna variable, menostat por ejemplo, la que mayor diferencia muestra a priori.

bc.fit.km.menostat<-survfit(surv.obj~menostat, data=BC, type="kaplan-meier")
head(fortify(bc.fit.km.menostat, strata=menostat))
##   time n.risk n.event n.censor      surv     std.err     upper     lower
## 1 0.03    621       1        0 0.9983897 0.001611604 1.0000000 0.9952411
## 2 0.16    620       2        0 0.9951691 0.002795893 1.0000000 0.9897306
## 3 0.30    618       1        0 0.9935588 0.003231035 0.9998707 0.9872867
## 4 0.33    617      13        0 0.9726248 0.006732249 0.9855436 0.9598753
## 5 0.53    604       1        0 0.9710145 0.006933169 0.9842994 0.9579089
## 6 0.62    603       1        0 0.9694042 0.007129068 0.9830445 0.9559532
##   strata
## 1   Post
## 2   Post
## 3   Post
## 4   Post
## 5   Post
## 6   Post
ggsurvplot(bc.fit.km,  xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability",  ggtheme = theme_bw())

ggsurvplot(bc.fit.km,  xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability",  ggtheme = theme_bw(), surv.median.line = "h", risk.table = T, cumevents=T)

Junto con la curva, también tenemos la opción de representar las tablas de riesgo y eventos de esta misma.

ggsurvplot(bc.fit.km.menostat, xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability", pval=TRUE,  ggtheme = theme_bw())

3.1.2.Desviación estandar e intervalos de confianza

Se definen directamente en la función survfit().

Se puede definir el error: tipo de estimación para las desviaciones, conf.type para definit el tipo de transformación para calcular los intervalos de confianza y conf.int, el nivel de confianza (95% por defecto) (Este nivel quiere decir que ese tanto porciento de las veces la media se encontrará dentro del intervalo).

  • error: greenwood(defecto), tsiatis

  • conf.type: plain, log(defecto), log-log

surv.tsiatis.plain<-survfit(surv.obj ~ 1, data=BC,  error="tsiatis", conf.type = "plain", conf.int = 0.99)
head(surv.tsiatis.plain)
## Call: survfit(formula = surv.obj ~ 1, data = BC, error = "tsiatis", 
##     conf.type = "plain", conf.int = 0.99)
## 
##       n  events  median 0.99LCL 0.99UCL 
##   817.0   738.0    27.8    23.7    32.3
surv.tsiatis.loglog<-survfit(surv.obj ~ 1, data=BC,  error="tsiatis", conf.type = "log-log", conf.int = 0.90)
head(surv.tsiatis.loglog)
## Call: survfit(formula = surv.obj ~ 1, data = BC, error = "tsiatis", 
##     conf.type = "log-log", conf.int = 0.9)
## 
##      n events median 0.9LCL 0.9UCL 
##  817.0  738.0   27.8   24.7   30.6

3.1.3.Media, mediana, percentiles

media = \(\int^{\tau}_{0}S(t)dt\) siendo S(t) la función de supervivencia estimada hasta valor \(\tau\).

print(bc.fit.km, print.rmean=T)
## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "kaplan-meier")
## 
##          n     events     *rmean *se(rmean)     median    0.95LCL 
##     817.00     738.00      40.62       1.27      27.83      24.47 
##    0.95UCL 
##      30.98 
##     * restricted mean with upper limit =  217
probs<-c(0.05, 0.5, 0.95)
quantile(bc.fit.km, probs)
## $quantile
##      5     50     95 
##   2.60  27.83 108.28 
## 
## $lower
##      5     50     95 
##   1.02  24.47 105.98 
## 
## $upper
##      5     50     95 
##   5.85  30.98 118.50
quantile(bc.fit.km.menostat, probs)
## $quantile
##                  5    50     95
## menostat=Post 1.77 24.11 107.85
## menostat=Pre  7.29 40.37 115.18
## 
## $lower
##                  5    50     95
## menostat=Post 0.99 22.08 101.54
## menostat=Pre  4.40 32.75 107.62
## 
## $upper
##                   5    50     95
## menostat=Post  5.32 27.56 118.36
## menostat=Pre  11.33 50.33     NA

3.1.4. Estimación de la función de riesgo acumulado

Nos permitirá ver como aumenta con el tiempo el riesgo de que ocurra el evento. Como es de esperar, veremos claramente como es inversa a la curva de supervivencia, ya que cuando el riesgo sube, esta baja.

Tenemos varias formas de estimarla, entre ellas exploraremos máxima verosimilitud y Nelson-Aalen

El riesgo acumulado vendrá dado por CumHaz= cumsum(n.event/n.risk), usando la función mutate (perteneciente a la libreria dplyr).Para visualizarlo simplemente tendremos que indicar fun="cumhaz".

Podemos visualizarlo para los diferentes niveles de un evento en concreto usando group_by.

3.1.4.1. Máxima verosimilitud

En este caso se estima igual que anteriormente, pero se calcula la suma acumulada del número de enventos entre las personas en riesgo.

bc.fit.ra<-survfit(surv.obj~1) %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
head(bc.fit.ra)
##   time n.risk n.event n.censor      surv     std.err     upper     lower
## 1 0.03    817       1        0 0.9987760 0.001224740 1.0000000 0.9963814
## 2 0.16    816       2        0 0.9963280 0.002123916 1.0000000 0.9921891
## 3 0.23    814       1        0 0.9951040 0.002453995 0.9999018 0.9903293
## 4 0.30    813       1        0 0.9938800 0.002745339 0.9992423 0.9885466
## 5 0.33    812      15        0 0.9755202 0.005542106 0.9861744 0.9649811
## 6 0.53    797       1        0 0.9742962 0.005682535 0.9852081 0.9635051
##        CumHaz
## 1 0.001223990
## 2 0.003674971
## 3 0.004903472
## 4 0.006133484
## 5 0.024606391
## 6 0.025861096
ggsurvplot(bc.fit.ra, fun="cumhaz", xlab="Tiempo(meses)", ylab="Riempo acumulado", main="Riesgo acumulado", censor=T, conf.int=T, palette="Dark2",  ggtheme = theme_bw())

bc.fit.ra.menostat<-survfit(surv.obj~menostat, data=BC) %>% fortify %>% group_by(strata) %>% mutate(CumHaz = cumsum(n.event/n.risk))
head(bc.fit.ra.menostat)
## # A tibble: 6 x 10
## # Groups:   strata [1]
##    time n.risk n.event n.censor  surv std.err upper lower strata  CumHaz
##   <dbl>  <dbl>   <dbl>    <dbl> <dbl>   <dbl> <dbl> <dbl> <fct>    <dbl>
## 1  0.03    621       1        0 0.998 0.00161 1     0.995 Post   0.00161
## 2  0.16    620       2        0 0.995 0.00280 1     0.990 Post   0.00484
## 3  0.3     618       1        0 0.994 0.00323 1.000 0.987 Post   0.00645
## 4  0.33    617      13        0 0.973 0.00673 0.986 0.960 Post   0.0275 
## 5  0.53    604       1        0 0.971 0.00693 0.984 0.958 Post   0.0292 
## 6  0.62    603       1        0 0.969 0.00713 0.983 0.956 Post   0.0308
ggsurvplot(survfit(surv.obj~menostat, data=BC), fun="cumhaz", xlab="Tiempo(meses)", ylab="Riesgo acumulado", main="Riesgo acumulado", censor=T, conf.int=T, legend.title="Perfil menostat", palette="Dark2",  ggtheme = theme_bw())

También podemos ver los valores de riesgo en una tabla:

ggsurvplot(bc.fit.km.menostat, xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability", pval=TRUE, risk.table = TRUE, risk.table.y.text.col = TRUE,  ggtheme = theme_bw())

3.1.4.2. Nelson-Aalaen

La información es extraída de la misma manera que la Estimación Máxima Verosimilitud, con la diferencia que se calcula la suma acumulada del número de eventos entre las personas en riesgo.

qplot(time, CumHaz, data=bc.fit.ra, geom="step", xlab="Tiempo(meses)", ylab="Riesgo acumulado", main="Riesgo acumulado")

qplot(time, CumHaz, col=strata, data=bc.fit.ra.menostat, geom="step", xlab="Tiempo(meses)", ylab="Riesgo acumulado", main="Riesgo acumulado")

3.2. Métodos no paramétricos

Todo lo visto anteriormente se realizada de forma no paramétrica, pero solo lo hemos visto para Kaplain-Meier.

Veamos si otros métodos nos dan alguna diferencia.

km<-survfit(surv.obj~1, data=BC, type="kaplan-meier")
fh<-survfit(surv.obj~1, data=BC, type="fleming-harrington")
fh2<-survfit(surv.obj~1, data=BC, type="fh2")
list.fit=list(km,fh)

ggplot() + geom_step(aes(time, surv, col = "Kaplain-meier"), data = km) + geom_step(aes(time, surv, col = "Fleming-harrington"), data = fh) + geom_step(aes(time, surv, col = "fh2"), data=fh2) + labs(x = "Tiempo (meses)", y = "Probabilidad de Supervivencia", col = "Ajustes no parametricos", title = "Comparacion de ajustes en curvas de supervivencia")

La diferencia de Kaplain-meier con otros ajustes es inapreciable a simple vista, por lo que seguiremos usando este primero.

3.3. Métodos paramétricos

Estas técnicas a diferencia de las no paramétricas, tienen que ajustarse a una distribución de datos, y cumplir condiciones de validez.

Para estimar las curvas de supervivencia de forma no paramétrica, usaremos la función flexsurvreg(), donde los parámetros se estiman por máxima verosimilitud (visto en 3.1.4.1).

Tenemos disponibles multitud de distribuciones por las que ajustar (indicadas en el parámetro dist), algunas de ellas son:

Nombre dist Numero de parámetros
Gamma generalizada (stable) “gengamma” 3
Gamma generalizada (original) “gengamma.orig” 3
F generalizada (stable) “genf” 4
F generalizada (original) “genf.orig” 4
Weibull “weibull” 2
Gamma “gamma” 2
Exponencial “exp” 1
Log-logistica “llogis” 2
Log-normal “lnorm” 2
Gompertz “gompertz” 2

Los parámetros definidos como positivos se estiman con la logarítmica y los intervalos de confianza se calculan a partir de la matriz Hessiana.

Compararemos graficamente algunas de ellas aplicadas a nuestro dataset:

#Distribuciones
dist<-c("exp", "gengamma", "gamma", "weibull", "llogis", "lnorm", "gompertz", "genf")
#Lista de ajustes de las distribuciones
fit.param<-sapply(dist, function(x) flexsurvreg(surv.obj~1, dist=x), USE.NAMES=T, simplify=F)

param.dist<-names(fit.param) %>% sapply(function(x) cbind(dist = x, data.frame(summary(fit.param[[x]]))), simplify = F) %>% do.call("rbind", .)

x <- (1:length(dist)) + 1  #Colores
names(x) <- dist  

#Supervivence
ggplot() + geom_line(aes(time, est, col = factor(dist)), param.dist, size=0.9, alpha=0.5) + labs(col = "Distribuciones parametricas", linetype = "Ajustes", x ="Tiempo(meses)", y = "Probabilidad de Supervivencia", title = "Comparacion de ajustes parametricos")

#Risk
ggplot() + geom_line(aes(time, -log(est), col = factor(dist)), param.dist, size=0.9, alpha=0.5) + labs(col = "Distribuciones parametricas", linetype = "Ajustes", x ="Tiempo(meses)", y = "Riesgo acumulado", title = "Comparacion de ajustes parametricos")

También podemos ver los valores de AIC, BIC y LogLik para cada distribución.

  • AIC (Akaike information criterion): Estimador de la calidad relativa del modelo estadístico para el dataset al que se ha aplicado. Representa la cantidad de información perdida, por tanto el modelo será mejor para un menor AIC.

  • BIC (Bayesian information criterion): Criterio para seleccionar alguno de los modelos, es preferible el de menos BIC.

  • LogLik (likelihood function)

AIC.model <- sapply(fit.param, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = logLik(x)), simplify = T)
AIC.model
##              exp  gengamma     gamma   weibull    llogis     lnorm
## AIC     6972.944  6957.689  6958.970  6956.217  7027.395  7155.689
## BIC     6977.649  6971.806  6968.382  6965.629  7036.806  7165.101
## LogLik -3485.472 -3475.844 -3477.485 -3476.109 -3511.697 -3575.845
##         gompertz      genf
## AIC     6959.958  6959.691
## BIC     6969.369  6978.513
## LogLik -3477.979 -3475.845

No hay demasiada diferencia entre los modelos, pero a simple vista el de weibull parace ser el más adecuado en este caso.

3.4.Regresión de Cox

El objetivo principal de Cox (modelo semiparamétrico de regresión multivariante) es poder estimar el efecto de riesgo conjunto de distintas variables, además, funciona con variables continuas además de categóricas, permitiendonos esto analizar otras variables no vistas anteriormente como por ejemplo la edad.

Podemos ver como cada factor influye al Hazard rate o failure rate. Como cambia la supervivencia respecto a la otra variable. Es la probabilidad de que si algo sobrevive a un momento, lo hará también para el siguiente.

Función de riesgo: \(h(t)=\frac{f(t)}{R(t)}\) (riesgo de que se produzca en evento en el instante de tiempo \(t\).)

Siendo:

  • \(f(t)\): Función de probabilidad de densidad (probabilidad de que el evento caiga en un intervalo específico).

  • \(R(t)\): Función de supervivencia dado por las covariantes \(x_1, x_2,... x_p\) indicadas, siendo estos multiplicados por sus respectivos impactos \(b_1, b_2,...b_p\).

\(h0\): Baseline hazard (no covariables, all \(x_i\)=0)

Hazard ratios (HR) \(exp(b_i)\)

  • HR=1: Esa variable no tiene efecto
  • HR <1: Reducción de la probabilidad hazard
  • HR >1: aumento

Prueba de hazards proporcionales o proporcionalidad de los riesgos:

“If an individual has a risk of death at some initial time point that is twice as high as that of another individual, then at all later times the risk of death remains twice as high.” Se realizarán con cox.zph para aceptar el modelo como válido, si esta da un \(p-value \ge 0.05\). De esta forma sabemos que se mantiene la diferencia entre las pacientes.

3.4.1. Busqueda exhaustiva

3.4.1.1. Modelos univariantes

library(FSelector)

returnSignificatives<-function(vbles) {
  variables<-list()
  pvalues<-list()
  for (v in vbles) {
    form<-as.simple.formula(v,"surv.obj")
    sum<-summary(coxph(form, data=BC))$coefficients
    p<-sum[,"Pr(>|z|)"]
    if(p<0.05){
      variables=c(variables,v)
      pvalues=c(pvalues,p)
    }
  }
  significancia<-data.frame(unlist(variables),unlist(pvalues))
  names(significancia)<-c("variable", "p-value")
  return(significancia)
}

returnSignificatives(c("age", "menostat","tgrade","pnodes", "Luminal", "NodalStatus", "NodeStatus"))
##   variable      p-value
## 1      age 7.483853e-05
## 2 menostat 1.362899e-04
## 3   pnodes 2.374573e-03

Parece ser que individualmente, las más significativas son edad, menopausia y número de nodos.

3.4.1.2. Modelos multivariantes

Probaremos juntas las más significativas, para ver si explican mejor la función de riesgo, pero algunas como NodeStatus o tgrade no las descartaremos aún ya que pueden dar buenos resultados.

Al añadir más variables, se tomará una de ellas como referencia.

Con la función de búsqueda exhaustiva intentaremos encontrar la mejor combinación posible de variables que expliquen la función de riesgo. Nos quedaremos en cada nivel con la mejor (la de menor \(p-value\)) e iremos añadiendo otras solo en el caso de que su incorporación mejore el resultado que ya teníamos. Nos quedamos con la evolución de los \(p-values\).

exhaustiveSearch <-function(lista.vbles) {
  
  pvalue.acumulado<-list()
  
  for(nivel in 1:length(lista.vbles)) { #cada nivel (#vbles=#nivel)
    best.pvalue<-1
    for(v in lista.vbles) { #vbles de cada nivel
      ifelse(nivel==1, att<-v, att<-c(v, lista.acumulada))
      ec <-as.simple.formula(att,"surv.obj")
      sum<-summary(coxph(ec, data=BC))
      new.p<-sum$logtest["pvalue"]
      if(new.p <0.05 && new.p<best.pvalue) {
        best.pvalue<-new.p
        best.vble.nivel<-v
      }
    }
    
    lista.vbles<-lista.vbles[lista.vbles != best.vble.nivel]
    
    ifelse(nivel!=1,  pvalue.acumulado<-c(pvalue.acumulado, best.pvalue), first.pvalue<-best.pvalue)
    ifelse(nivel==1, lista.acumulada<-best.vble.nivel,lista.acumulada<-c(lista.acumulada, best.vble.nivel))
    ifelse(nivel==1, p.acum.ac<-best.pvalue, p.acum.ac<-pvalue.acumulado[nivel-1])
    if(length(lista.acumulada)!=nivel || p.acum.ac>best.pvalue[1] || length(lista.acumulada)==10) break
  }
  
df<-data.frame(lista.acumulada,unlist(list(c(first.pvalue,pvalue.acumulado))))
names(df)<-c("vble acumulada", "pvalue_acumulado")
return(df)
}

es<-exhaustiveSearch(c("age", "menostat","tgrade","pnodes", "Luminal", "NodalStatus", "NodalStatus2", "NodeStatus", "GrupoEdadCluster", "GrupoEdad2"))

df.order<-es[order(es$pvalue_acumulado),]
df.order
##      vble acumulada pvalue_acumulado
## 3            tgrade     5.956164e-09
## 4            pnodes     9.276238e-09
## 5        GrupoEdad2     1.690446e-08
## 2          menostat     3.034877e-08
## 6               age     3.851068e-08
## 7  GrupoEdadCluster     8.173465e-08
## 8        NodeStatus     1.851588e-07
## 9       NodalStatus     3.260355e-07
## 10          Luminal     2.267670e-06
## 1      NodalStatus2     6.423682e-06

Las variables más influyentes son tgrade, pnodes y GrupoEdad. Efectivamente coincide con las variables encontrados en el análisis multivariante, pero ahora el pvalue es mucho menor, pasando de ser del orden \(e-05, e-04, e-03\) a \(e-09\), lo cual nos confirma que juntas explican mejor la función de rieso.

En el primer puesto tenemos a tgrade, variable que no aparecia anteriormente. Sola no es significativa, pero sí en conjunto.


Así como para añadir una variable a analizar se usa +, para ver todas las interacciones posibles entre las variables usaremos *.

Comprobamos que sea válido antes de examinarlo.

El coxph() nos dará la siguiente información:

  • coef: Coeficiente estimado de la beta (los valores beta, en función del valor explican la importancia de esta variable para la supervivencia)
  • exp(coef): Exponencial elevado al coeficiente beta estimado.
  • se(coef): La desviación estándar de la estimación.
  • z: Valor del estadístico para prueba Wald de beta igual a cero.
  • p: p-valor de la prueba Wald, beta igual a cero.

Usando summary() los resultados son una versión más detallada de la función coxph, en esta se encuentran agregada una tabla con los coeficientes estimados de los exponentes β con sus respectivos intervalos de confianza.

Además se agregan estadísticos de Concordancia y R2; así como los estadísticos de prueba, grados de libertad y p-valores de las pruebas de hipótesis de Razón de Verosimilitud, Prueba de Wald y Prueba de Puntajes (score test).

res.cox <- coxph(surv.obj~tgrade+pnodes+GrupoEdadCluster, data=BC)
cox.zph(res.cox, transform="km", global=TRUE) #ACEPTADO
##                              rho   chisq      p
## tgradeII                -0.10345  7.8452 0.0051
## tgradeIII               -0.06237  2.9899 0.0838
## tgradeIV                -0.05267  2.0735 0.1499
## pnodes                   0.00662  0.0288 0.8653
## GrupoEdadCluster(50,66] -0.03672  0.9831 0.3214
## GrupoEdadCluster(66,90] -0.07404  4.0399 0.0444
## GLOBAL                        NA 12.1104 0.0596
summ<-summary(res.cox)
summ
## Call:
## coxph(formula = surv.obj ~ tgrade + pnodes + GrupoEdadCluster, 
##     data = BC)
## 
##   n= 817, number of events= 738 
## 
##                             coef exp(coef) se(coef)      z Pr(>|z|)    
## tgradeII                 0.09981   1.10496  0.09961  1.002 0.316336    
## tgradeIII                0.05582   1.05741  0.14269  0.391 0.695658    
## tgradeIV                -0.98886   0.37200  0.51252 -1.929 0.053678 .  
## pnodes                   0.02798   1.02837  0.01089  2.570 0.010182 *  
## GrupoEdadCluster(50,66]  0.28066   1.32401  0.08752  3.207 0.001342 ** 
## GrupoEdadCluster(66,90]  0.37929   1.46125  0.10241  3.703 0.000213 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## tgradeII                    1.105     0.9050    0.9090     1.343
## tgradeIII                   1.057     0.9457    0.7994     1.399
## tgradeIV                    0.372     2.6882    0.1362     1.016
## pnodes                      1.028     0.9724    1.0067     1.051
## GrupoEdadCluster(50,66]     1.324     0.7553    1.1153     1.572
## GrupoEdadCluster(66,90]     1.461     0.6843    1.1955     1.786
## 
## Concordance= 0.578  (se = 0.012 )
## Rsquare= 0.037   (max possible= 1 )
## Likelihood ratio test= 31.06  on 6 df,   p=2e-05
## Wald test            = 30.73  on 6 df,   p=3e-05
## Score (logrank) test = 31.13  on 6 df,   p=2e-05
res.cox1 <- coxph(surv.obj~tgrade+NodeStatus+NodalStatus, data=BC)
  cox.zph(res.cox1, transform="km", global=TRUE) #DESCARTADO
##                               rho    chisq      p
## tgradeII                -0.091220 6.01e+00 0.0142
## tgradeIII               -0.025336 5.07e-01 0.4767
## tgradeIV                -0.039975 1.22e+00 0.2702
## NodeStatusNode positive  0.020980 3.45e-01 0.5571
## NodalStatusG1           -0.023330 4.32e-01 0.5112
## NodalStatusG2           -0.000196 2.98e-05 0.9956
## NodalStatusG3           -0.084237 5.23e+00 0.0222
## NodalStatusG4                  NA      NaN    NaN
## GLOBAL                         NA 2.35e+01 0.0028
res.cox2 <- coxph(surv.obj~menostat+tgrade+NodeStatus+Luminal, data=BC) 
cox.zph(res.cox2, transform="km", global=TRUE) #DESCARTADO
##                               rho    chisq      p
## menostatPre              0.094968 6.58e+00 0.0103
## tgradeII                -0.080805 4.77e+00 0.0289
## tgradeIII               -0.043925 1.43e+00 0.2314
## tgradeIV                -0.049536 1.84e+00 0.1755
## NodeStatusNode positive -0.045546 1.55e+00 0.2130
## LuminalLuminalA          0.010433 8.15e-02 0.7753
## LuminalLuminalB          0.026295 5.13e-01 0.4737
## LuminalTN               -0.000543 2.19e-04 0.9882
## GLOBAL                         NA 1.70e+01 0.0300
res.cox3 <- coxph(surv.obj~menostat+tgrade+NodeStatus+Luminal+GrupoEdadCluster, data=BC) 
cox.zph(res.cox3, transform="km", global=TRUE) #ACEPTADO
##                              rho   chisq      p
## menostatPre              0.07173  3.4189 0.0645
## tgradeII                -0.08422  5.2189 0.0223
## tgradeIII               -0.04750  1.6938 0.1931
## tgradeIV                -0.05252  2.0714 0.1501
## NodeStatusNode positive -0.04050  1.2398 0.2655
## LuminalLuminalA          0.00593  0.0267 0.8701
## LuminalLuminalB          0.02757  0.5663 0.4517
## LuminalTN               -0.00431  0.0139 0.9062
## GrupoEdadCluster(50,66]  0.01945  0.2555 0.6132
## GrupoEdadCluster(66,90] -0.01017  0.0709 0.7900
## GLOBAL                        NA 17.9113 0.0565
summ3<-summary(res.cox3)
summ3
## Call:
## coxph(formula = surv.obj ~ menostat + tgrade + NodeStatus + Luminal + 
##     GrupoEdadCluster, data = BC)
## 
##   n= 817, number of events= 738 
## 
##                              coef exp(coef)  se(coef)      z Pr(>|z|)  
## menostatPre             -0.213294  0.807918  0.122083 -1.747   0.0806 .
## tgradeII                 0.059617  1.061430  0.106317  0.561   0.5750  
## tgradeIII                0.137181  1.147036  0.141452  0.970   0.3321  
## tgradeIV                -0.959077  0.383247  0.513807 -1.867   0.0620 .
## NodeStatusNode positive  0.125685  1.133925  0.089555  1.403   0.1605  
## LuminalLuminalA         -0.149212  0.861387  0.193617 -0.771   0.4409  
## LuminalLuminalB         -0.004047  0.995961  0.212315 -0.019   0.9848  
## LuminalTN               -0.120004  0.886917  0.209885 -0.572   0.5675  
## GrupoEdadCluster(50,66]  0.171218  1.186750  0.118177  1.449   0.1474  
## GrupoEdadCluster(66,90]  0.246775  1.279892  0.134044  1.841   0.0656 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## menostatPre                0.8079     1.2377    0.6360     1.026
## tgradeII                   1.0614     0.9421    0.8618     1.307
## tgradeIII                  1.1470     0.8718    0.8693     1.513
## tgradeIV                   0.3832     2.6093    0.1400     1.049
## NodeStatusNode positive    1.1339     0.8819    0.9514     1.351
## LuminalLuminalA            0.8614     1.1609    0.5894     1.259
## LuminalLuminalB            0.9960     1.0041    0.6569     1.510
## LuminalTN                  0.8869     1.1275    0.5878     1.338
## GrupoEdadCluster(50,66]    1.1867     0.8426    0.9414     1.496
## GrupoEdadCluster(66,90]    1.2799     0.7813    0.9842     1.664
## 
## Concordance= 0.577  (se = 0.012 )
## Rsquare= 0.039   (max possible= 1 )
## Likelihood ratio test= 32.64  on 10 df,   p=3e-04
## Wald test            = 30.93  on 10 df,   p=6e-04
## Score (logrank) test = 31.32  on 10 df,   p=5e-04
#cox1
summ$logtest[3]
##       pvalue 
## 2.468529e-05
#cox3
summ3$logtest[3]
##       pvalue 
## 0.0003134935
res.cox4 <- coxph(surv.obj~tgrade+NodeStatus+Luminal+GrupoEdadCluster, data=BC) 
cox.zph(res.cox4, transform="km", global=TRUE) #ACEPTADO
##                              rho   chisq      p
## tgradeII                -0.08405  5.1935 0.0227
## tgradeIII               -0.04689  1.6479 0.1992
## tgradeIV                -0.05256  2.0742 0.1498
## NodeStatusNode positive -0.03988  1.2030 0.2727
## LuminalLuminalA          0.00814  0.0505 0.8223
## LuminalLuminalB          0.02946  0.6455 0.4217
## LuminalTN               -0.00422  0.0133 0.9082
## GrupoEdadCluster(50,66] -0.03870  1.1224 0.2894
## GrupoEdadCluster(66,90] -0.07259  3.8113 0.0509
## GLOBAL                        NA 15.0513 0.0895
summ4<-summary(res.cox4)
summ4$logtest[3]
##       pvalue 
## 0.0005183621
Variables p-value propircionalidad p-value modelo
tgrade+pnodes+GrupoEdadCluster 0.0596 2.468529e-05
menostat+tgrade+NodeStatus+Luminal+GrupoEdadCluster 0.0565 0.0003134935
tgrade+NodeStatus+Luminal+GrupoEdadCluster 0.0895 0.0005183621

Algunos daban un buen p-value como modelo pero hemos tenido que descartarlo ya que no cumplia la proporcionalidad. De entre los restantes probados, el segundo es el más fiable en cuanto a proporcionalidad, pero con poca diferencia respecto al primero, el cual marca una notable diferencia en cuanto al p-value del modelo de cox.

El resultado de proporcionalidad mejora mucho al quitar menostat, pero no el pvalue del modelo.

Probaremos alguno más siguiendo el patrón con mejores resultados.

res.cox5 <- coxph(surv.obj~menostat+tgrade+NodeStatus+NodalStatus+Luminal+GrupoEdadCluster, data=BC) 
cox.zph(res.cox5, transform="km", global=TRUE) #RECHAZADO
##                              rho    chisq       p
## menostatPre              0.06805  2.97541 0.08454
## tgradeII                -0.09666  6.85358 0.00885
## tgradeIII               -0.03388  0.93281 0.33413
## tgradeIV                -0.04972  1.89922 0.16817
## NodeStatusNode positive  0.02557  0.52991 0.46664
## NodalStatusG1           -0.03067  0.76689 0.38118
## NodalStatusG2           -0.00489  0.01898 0.89041
## NodalStatusG3           -0.08658  5.60686 0.01789
## NodalStatusG4                 NA      NaN     NaN
## LuminalLuminalA          0.00139  0.00147 0.96940
## LuminalLuminalB          0.03365  0.84761 0.35723
## LuminalTN               -0.00643  0.03108 0.86007
## GrupoEdadCluster(50,66]  0.01962  0.25098 0.61639
## GrupoEdadCluster(66,90] -0.01056  0.07693 0.78150
## GLOBAL                        NA 31.41162 0.00485

Al añador el grupo nodal se anula el criterio de proporcionalidad.

res.cox6 <- coxph(surv.obj~menostat+NodeStatus+Luminal+GrupoEdadCluster, data=BC) 
cox.zph(res.cox6, transform="km", global=TRUE) #ACEPTADO
##                              rho    chisq      p
## menostatPre              0.07539  3.78033 0.0519
## NodeStatusNode positive -0.07990  4.74086 0.0295
## LuminalLuminalA          0.01216  0.11216 0.7377
## LuminalLuminalB          0.02473  0.45645 0.4993
## LuminalTN               -0.00315  0.00743 0.9313
## GrupoEdadCluster(50,66]  0.01726  0.20065 0.6542
## GrupoEdadCluster(66,90] -0.00680  0.03171 0.8587
## GLOBAL                        NA 12.09198 0.0976
summ6<-summary(res.cox6)

Al quitar el grado mejora mucho la proporcionalidad.

survfit(surv.obj ~ tgrade, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grado", conf.int = T, legend.title = "Grado", pval=TRUE,  ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.32)

Si visualizamos su curva de supervivencia podemos ver como tiene un pvalue poco fiable.El nodo IV tiene la menor supervivencia como podemos ver muy claramente, pero al haber tan pocos pacientes en ese grupo sus intervalos de confianza son muy amplios para poder alcanzar el 95% que le hemos marcado por defecto, solapando completamente al resto de niveles.

res.cox7 <- coxph(surv.obj~menostat+NodeStatus+NodalStatus+Luminal+GrupoEdadCluster, data=BC) 
cox.zph(res.cox7, transform="km", global=TRUE) #ACEPTADO
## Warning in cor(xx, r2): the standard deviation is zero
##                              rho   chisq      p
## menostatPre              0.07219  3.4030 0.0651
## NodeStatusNode positive  0.02766  0.5723 0.4493
## NodalStatusG1           -0.04701  1.6636 0.1971
## NodalStatusG2           -0.01631  0.1994 0.6552
## NodalStatusG3           -0.09283  6.4465 0.0111
## NodalStatusG4                 NA     NaN    NaN
## LuminalLuminalA          0.00649  0.0321 0.8578
## LuminalLuminalB          0.02755  0.5680 0.4511
## LuminalTN               -0.00547  0.0225 0.8809
## GrupoEdadCluster(50,66]  0.01659  0.1826 0.6691
## GrupoEdadCluster(66,90] -0.00624  0.0268 0.8699
## GLOBAL                        NA 23.9073 0.0131
summ7<-summary(res.cox7)
res.cox8 <- coxph(surv.obj~NodeStatus+NodalStatus+Luminal, data=BC) 
summ8<-summary(res.cox8)

Habiendo quitado el grado tumoral que causaba problemas, se puede añadir NodalStatus y la propocionalidad mejora bastante.

summ4$logtest[3]
##       pvalue 
## 0.0005183621
summ6$logtest[3]
##       pvalue 
## 0.0006044408
summ7$logtest[3]
##       pvalue 
## 0.0002492881
summ8$logtest[3]
##    pvalue 
## 0.0430766

A diferencia de los demás que no marcan gran diferencia, quitarle la información sobre el grado tumoral el resultado mejora enormemente, además, la edad y estado menopausico mejoran mucho el pvalue, por lo que parecen ser importantes en el estudio. Al añadirle otros empeora un poco el resultado por lo que en principio menostat+NodeStatus+NodalStatus+Luminal+GrupoEdadCluster serán neustras covariantes principales.

Exploramos varios valores y finalmente podemos afinar los levels que mejor explican la función de riesgo, dando el menor p-value encontrado ahasta el momento:

res.cox8<-coxph(surv.obj~+I(menostat=="Post")+ I(NodeStatus=="NodePositive") + I(NodalStatus=="G2")+I(GrupoEdadCluster=="(50,66]")+I(Luminal=="TN"), data=BC)
summ8<-summary(res.cox8)
summ8$logtest[3]
##      pvalue 
## 0.003594444
summ7
## Call:
## coxph(formula = surv.obj ~ menostat + NodeStatus + NodalStatus + 
##     Luminal + GrupoEdadCluster, data = BC)
## 
##   n= 817, number of events= 738 
## 
##                             coef exp(coef) se(coef)      z Pr(>|z|)  
## menostatPre             -0.22285   0.80023  0.12101 -1.842   0.0655 .
## NodeStatusNode positive  0.31993   1.37703  0.18799  1.702   0.0888 .
## NodalStatusG1           -0.31326   0.73106  0.20566 -1.523   0.1277  
## NodalStatusG2           -0.26269   0.76898  0.20506 -1.281   0.2002  
## NodalStatusG3           -0.01144   0.98863  0.19368 -0.059   0.9529  
## NodalStatusG4                 NA        NA  0.00000     NA       NA  
## LuminalLuminalA         -0.13895   0.87027  0.19366 -0.717   0.4731  
## LuminalLuminalB         -0.04109   0.95974  0.21239 -0.193   0.8466  
## LuminalTN               -0.10574   0.89966  0.21011 -0.503   0.6148  
## GrupoEdadCluster(50,66]  0.13973   1.14996  0.11722  1.192   0.2332  
## GrupoEdadCluster(66,90]  0.19479   1.21505  0.13409  1.453   0.1463  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## menostatPre                0.8002     1.2496    0.6313     1.014
## NodeStatusNode positive    1.3770     0.7262    0.9526     1.990
## NodalStatusG1              0.7311     1.3679    0.4885     1.094
## NodalStatusG2              0.7690     1.3004    0.5145     1.149
## NodalStatusG3              0.9886     1.0115    0.6763     1.445
## NodalStatusG4                  NA         NA        NA        NA
## LuminalLuminalA            0.8703     1.1491    0.5954     1.272
## LuminalLuminalB            0.9597     1.0419    0.6329     1.455
## LuminalTN                  0.8997     1.1115    0.5960     1.358
## GrupoEdadCluster(50,66]    1.1500     0.8696    0.9139     1.447
## GrupoEdadCluster(66,90]    1.2151     0.8230    0.9342     1.580
## 
## Concordance= 0.584  (se = 0.012 )
## Rsquare= 0.04   (max possible= 1 )
## Likelihood ratio test= 33.23  on 10 df,   p=2e-04
## Wald test            = 33.3  on 10 df,   p=2e-04
## Score (logrank) test = 33.59  on 10 df,   p=2e-04

Los Luminals y Grupos nodales tienen el valor ß negativo , lo cual quiere decir que a mayor sea este, menor es el riesgo, los grupos de edad por el contrario tiene un ß positivo, a mayor número de nodos mayor riesgo, y de esta forma interpretamos el riesgo.

Los evaluaremos más adelante, pero podemos intuir que el riesgo principalmente para el estado nodal positivo.

También volvemos a ver que la (pre)menopausia influye enormemente, llegando a ser la de menor p-value.

El intervalo de confianza con mayor rango lo tenemos en el estado nodal.

Otro efecto que también hemos observado es que a mayor edad, menor es el riesgo.

3.4.2. Comparación modelos

La comparación entre los modelos se realiza con la función anova(), quw realiza la prueba de Razón de Verosimilitud entre dos modelo, le daremos como argumentos de la función son los dos modelos ajustados a comparar.

anova(res.cox8, res.cox7)
## Analysis of Deviance Table
##  Cox model: response is  surv.obj
##  Model 1: ~ +I(menostat == "Post") + I(NodeStatus == "NodePositive") + I(NodalStatus == "G2") + I(GrupoEdadCluster == "(50,66]") + I(Luminal == "TN")
##  Model 2: ~ menostat + NodeStatus + NodalStatus + Luminal + GrupoEdadCluster
##    loglik  Chisq Df P(>|Chi|)   
## 1 -4223.7                       
## 2 -4214.8 17.622  6   0.00725 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(res.cox7, res.cox6)
## Analysis of Deviance Table
##  Cox model: response is  surv.obj
##  Model 1: ~ menostat + NodeStatus + NodalStatus + Luminal + GrupoEdadCluster
##  Model 2: ~ menostat + NodeStatus + Luminal + GrupoEdadCluster
##    loglik  Chisq Df P(>|Chi|)  
## 1 -4214.8                      
## 2 -4218.7 7.6729  3   0.05328 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En la tabla se muestra para los dos modelos la log-verosimilitud, el estadísticos Chisq, los grados de libertad Df y el p-valor P(>|Chi|) de la prueba.


3.5. Análisis exploratorio

3.5.1. Distrbución covariables con más significancia

Distribución de los eventos en el tiempo:

ggsurvevents(fit=survfit(surv.obj ~ 1, BC), normalized=T)

Las pacientes en las que tiene ocurrencia el evento son la inmensa mayoria, distribuidos estos uniformemente en el tiempo, con una diferencia de altura no demasiado grande para unos 25 meses.

Además, teniendo ya una previa idea de la importancia afecta cada covariable, haremos un análisis y lo veremos gráficamente.

Lo haremos de forma no paramétrica con Kaplain-Meier, ya que los demás no marcan diferencia respecto a este.

Empezaremos viendo más concretamente que antes la distribución de los datos marcados como más significativos (En %):

with(BC, prop.table(table(menostat, cens)))
##         cens
## menostat          0          1
##     Post 0.08200734 0.67809058
##     Pre  0.01468788 0.22521420
with(BC, prop.table(table(NodalStatus, cens)))
##            cens
## NodalStatus          0          1
##          G0 0.02692778 0.40881273
##          G1 0.01468788 0.12239902
##          G2 0.01223990 0.12239902
##          G3 0.03182375 0.21052632
##          G4 0.01101591 0.03916769
with(BC, prop.table(table(Luminal, cens)))
##                cens
## Luminal                   0           1
##   HER2-enriched 0.006119951 0.035495716
##   LuminalA      0.055079559 0.603427173
##   LuminalB      0.014687882 0.119951040
##   TN            0.020807834 0.144430845
with(BC, table(Luminal, cens))
##                cens
## Luminal           0   1
##   HER2-enriched   5  29
##   LuminalA       45 493
##   LuminalB       12  98
##   TN             17 118
with(BC, prop.table(table(NodeStatus, cens)))
##                cens
## NodeStatus               0          1
##   Node negative 0.02692778 0.40881273
##   Node positive 0.06976744 0.49449204
with(BC, table(NodeStatus, cens))
##                cens
## NodeStatus        0   1
##   Node negative  22 334
##   Node positive  57 404

En principio podemos ver que de las muertes acontecidas, las pacientes post menopausicas ocupan la gran mayoría, a su vez, en este mismo grupo los eventos no observados son muy pocos (un 0,008%).

También comprobamos lo que hemos visto en el modelo de Cox, además de Grado 0, las pacientes con G3 son aquellas en las que más muertes se han observado, y G4 en las que menos. De hecho si volvemos a los resultados de Cox, este nivel lo encontrabamos con NA.

La distribución de pacientes para estado nodel es más o menos homogenea, viendo unas pocas más de muertes en estado nodal positivo, aunque también son más las pacientes en esta categoría.

Solo con estas tablas no podemos decir demasiado, veámoslo con otro tipo de análisis.

Luminal tiene casi todas sus pacientes concentradas en LuminalA. El mayor número de muertes está en LuminalA y Triple negativo.

3.5.2. Análisis exploratorio no-paramétrico

survfit(surv.obj ~ menostat, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por menopausia", conf.int = T, legend.title = "Menopausia", pval=TRUE,  ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.25)

survfit(surv.obj ~ NodeStatus, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por estado nodal", conf.int = T, legend.title = "Estado Nodal", pval=TRUE,  ggtheme = theme_bw(),  palette = c("#E7B800", "#2E9FDF"),  risk.table=T, risk.table.col="strata", risk.table.height = 0.25)

survfit(surv.obj ~ Luminal, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grado", conf.int = T, legend.title = "Grado", pval=TRUE,  ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.32)

survfit(surv.obj ~ NodalStatus, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grupo de nodos", conf.int = T, legend.title = "Estado de los nodos", pval=TRUE,  ggtheme = theme_bw(), palette = "Dark2", risk.table=T, risk.table.col="strata", risk.table.height = 0.35)

Otra vez más vemos que el p-value de las diferencias entre menopausia marca la mayor diferencia entre esos grupos con notable diferencia al resto de variables, teniendo las pre mejor supervivencia.

En estado nodal la supervivencia para Nodo positivo es algo menor, pero no es tanta la diferencia, además, sus intervalos de confianza se solapan.

Con un pvalue también muy bajo pero no tanto como menopausia, encontramos el grupo del estado nodal. Aun con una buena diferencia según el pvalue, los intervalos de confianza son grandes y no podemos apreciar mucho, tal vez decir que al principio el grupo G4 y G3 tienen una supervivencia menor (como era de esperar) y los grupos G0, G1 y G2 mayor.

Para todas coincide que a partir del tiempo 100 quedan muy pocas pacientes vivas.

Con el menos fiable de todos los pvalues tenemos Luminal. Como habiamos visto. HER2-Herinched al tener tan pocos pacientes en ese grupo sus intervalos de confianza son muy amplios para poder alcanzar el 95% que le hemos marcado por defecto.

Para seguir con el estudio eliminaremos este grupo que causa problemas y dificultades.

BC<-BC[as.factor(BC$Luminal)!="HER2-enriched",]
BC$Luminal<-factor(BC$Luminal, labels=c("LuminalA", "LuminalB", "TN" ))
surv.obj<-Surv(BC$time, BC$cens) 
survfit(surv.obj ~ Luminal, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grado", conf.int = T, legend.title = "Grado", pval=TRUE,  ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.32)

El \(pvalue\) mejora, aunque no hemos conseguido que los intervalos dejen de estar solapados. LuminalA tiene la mayor parte de los valores, pero B y TN al tener menos sus intervalos son mayores. Por lo visto graficamente parece que la peor supervivencia es para LuminalB, y a contrario de lo que habiamos previsto, la mejor para TN. No son unas deduciones muy fiables ya que como vemos en la tabla de riesgo, que queden mmuchas menos pacientes vivas e debido a que en un principio habia menos que en LuminalA.

Las colas largas independientes que vemos en todos los gráficos, como bien se puede ver en las tabals de riesgo son debidas a que solo llega un paciente vivo a los últimos intervalos medidos. En el caso de grupo nodal solo sobrevive G2, para LuminalA, Nodalmente positivo y Post menopausica.

En cuanto a los grupos de edad:

survfit(surv.obj ~ GrupoEdad2, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Grupos de edad", conf.int = T, legend.title = "Grupo edad", pval=TRUE,  ggtheme = theme_bw())

survfit(surv.obj ~ GrupoEdadCluster, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Grupos de edad Cluster", conf.int = T, legend.title = "Grupo edad Clustererizados", pval=TRUE,  ggtheme = theme_bw())

Así como tgrade dejaremos de explorarlo que que demuestra no ser una aportación importante en el estudio, vemos que los grupos de edad si marcan mucha diferencia, siendo mejor y más fiables los clusterizados ya que sus intervalos son menores los agrupados que los hechos arbitrariamente. Los usaremos para seguir viendo diferencias.

Si aplicamos el modelo de cox para ver la significancia de todas las interacciones posibles:

summary(coxph(surv.obj~menostat*NodalStatus*Luminal, data=BC)) 
## Call:
## coxph(formula = surv.obj ~ menostat * NodalStatus * Luminal, 
##     data = BC)
## 
##   n= 783, number of events= 709 
## 
##                                               coef exp(coef) se(coef)
## menostatPre                               -0.27061   0.76291  0.17232
## NodalStatusG1                              0.05886   1.06063  0.16337
## NodalStatusG2                              0.09613   1.10090  0.15637
## NodalStatusG3                              0.60381   1.82908  0.13382
## NodalStatusG4                              0.21456   1.23931  0.34247
## LuminalLuminalB                            0.44499   1.56047  0.20105
## LuminalTN                                  0.12469   1.13279  0.15610
## menostatPre:NodalStatusG1                  0.14013   1.15043  0.31010
## menostatPre:NodalStatusG2                  0.08851   1.09254  0.31403
## menostatPre:NodalStatusG3                 -0.57836   0.56082  0.28163
## menostatPre:NodalStatusG4                  0.59283   1.80911  0.46695
## menostatPre:LuminalLuminalB               -0.38681   0.67922  0.48097
## menostatPre:LuminalTN                     -0.02697   0.97339  0.30969
## NodalStatusG1:LuminalLuminalB             -0.45390   0.63515  0.41551
## NodalStatusG2:LuminalLuminalB             -0.38070   0.68338  0.41260
## NodalStatusG3:LuminalLuminalB             -0.86291   0.42193  0.30981
## NodalStatusG4:LuminalLuminalB             -0.15493   0.85647  0.59333
## NodalStatusG1:LuminalTN                   -0.30393   0.73791  0.46074
## NodalStatusG2:LuminalTN                    0.16604   1.18062  0.41010
## NodalStatusG3:LuminalTN                   -0.18090   0.83452  0.30305
## NodalStatusG4:LuminalTN                   -0.68125   0.50598  0.79733
## menostatPre:NodalStatusG1:LuminalLuminalB  0.79531   2.21512  0.74307
## menostatPre:NodalStatusG2:LuminalLuminalB -0.07554   0.92724  0.95515
## menostatPre:NodalStatusG3:LuminalLuminalB  1.28083   3.59962  0.67029
## menostatPre:NodalStatusG4:LuminalLuminalB       NA        NA  0.00000
## menostatPre:NodalStatusG1:LuminalTN       -0.47980   0.61891  0.76252
## menostatPre:NodalStatusG2:LuminalTN       -1.01026   0.36412  0.88979
## menostatPre:NodalStatusG3:LuminalTN        0.35075   1.42013  0.83802
## menostatPre:NodalStatusG4:LuminalTN        3.09207  22.02271  1.34340
##                                                z Pr(>|z|)    
## menostatPre                               -1.570  0.11633    
## NodalStatusG1                              0.360  0.71861    
## NodalStatusG2                              0.615  0.53870    
## NodalStatusG3                              4.512 6.42e-06 ***
## NodalStatusG4                              0.626  0.53099    
## LuminalLuminalB                            2.213  0.02687 *  
## LuminalTN                                  0.799  0.42443    
## menostatPre:NodalStatusG1                  0.452  0.65135    
## menostatPre:NodalStatusG2                  0.282  0.77807    
## menostatPre:NodalStatusG3                 -2.054  0.04001 *  
## menostatPre:NodalStatusG4                  1.270  0.20423    
## menostatPre:LuminalLuminalB               -0.804  0.42127    
## menostatPre:LuminalTN                     -0.087  0.93061    
## NodalStatusG1:LuminalLuminalB             -1.092  0.27467    
## NodalStatusG2:LuminalLuminalB             -0.923  0.35617    
## NodalStatusG3:LuminalLuminalB             -2.785  0.00535 ** 
## NodalStatusG4:LuminalLuminalB             -0.261  0.79400    
## NodalStatusG1:LuminalTN                   -0.660  0.50948    
## NodalStatusG2:LuminalTN                    0.405  0.68557    
## NodalStatusG3:LuminalTN                   -0.597  0.55056    
## NodalStatusG4:LuminalTN                   -0.854  0.39288    
## menostatPre:NodalStatusG1:LuminalLuminalB  1.070  0.28449    
## menostatPre:NodalStatusG2:LuminalLuminalB -0.079  0.93696    
## menostatPre:NodalStatusG3:LuminalLuminalB  1.911  0.05602 .  
## menostatPre:NodalStatusG4:LuminalLuminalB     NA       NA    
## menostatPre:NodalStatusG1:LuminalTN       -0.629  0.52920    
## menostatPre:NodalStatusG2:LuminalTN       -1.135  0.25621    
## menostatPre:NodalStatusG3:LuminalTN        0.419  0.67555    
## menostatPre:NodalStatusG4:LuminalTN        2.302  0.02135 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                           exp(coef) exp(-coef) lower .95
## menostatPre                                  0.7629    1.31076   0.54424
## NodalStatusG1                                1.0606    0.94283   0.77003
## NodalStatusG2                                1.1009    0.90834   0.81031
## NodalStatusG3                                1.8291    0.54672   1.40709
## NodalStatusG4                                1.2393    0.80690   0.63339
## LuminalLuminalB                              1.5605    0.64083   1.05226
## LuminalTN                                    1.1328    0.88277   0.83421
## menostatPre:NodalStatusG1                    1.1504    0.86924   0.62647
## menostatPre:NodalStatusG2                    1.0925    0.91530   0.59038
## menostatPre:NodalStatusG3                    0.5608    1.78311   0.32292
## menostatPre:NodalStatusG4                    1.8091    0.55276   0.72443
## menostatPre:LuminalLuminalB                  0.6792    1.47227   0.26461
## menostatPre:LuminalTN                        0.9734    1.02733   0.53050
## NodalStatusG1:LuminalLuminalB                0.6351    1.57444   0.28131
## NodalStatusG2:LuminalLuminalB                0.6834    1.46331   0.30441
## NodalStatusG3:LuminalLuminalB                0.4219    2.37004   0.22990
## NodalStatusG4:LuminalLuminalB                0.8565    1.16758   0.26771
## NodalStatusG1:LuminalTN                      0.7379    1.35518   0.29910
## NodalStatusG2:LuminalTN                      1.1806    0.84701   0.52848
## NodalStatusG3:LuminalTN                      0.8345    1.19829   0.46077
## NodalStatusG4:LuminalTN                      0.5060    1.97635   0.10603
## menostatPre:NodalStatusG1:LuminalLuminalB    2.2151    0.45144   0.51629
## menostatPre:NodalStatusG2:LuminalLuminalB    0.9272    1.07847   0.14262
## menostatPre:NodalStatusG3:LuminalLuminalB    3.5996    0.27781   0.96763
## menostatPre:NodalStatusG4:LuminalLuminalB        NA         NA        NA
## menostatPre:NodalStatusG1:LuminalTN          0.6189    1.61576   0.13886
## menostatPre:NodalStatusG2:LuminalTN          0.3641    2.74631   0.06366
## menostatPre:NodalStatusG3:LuminalTN          1.4201    0.70416   0.27479
## menostatPre:NodalStatusG4:LuminalTN         22.0227    0.04541   1.58257
##                                           upper .95
## menostatPre                                  1.0694
## NodalStatusG1                                1.4609
## NodalStatusG2                                1.4957
## NodalStatusG3                                2.3776
## NodalStatusG4                                2.4249
## LuminalLuminalB                              2.3141
## LuminalTN                                    1.5382
## menostatPre:NodalStatusG1                    2.1126
## menostatPre:NodalStatusG2                    2.0218
## menostatPre:NodalStatusG3                    0.9740
## menostatPre:NodalStatusG4                    4.5179
## menostatPre:LuminalLuminalB                  1.7435
## menostatPre:LuminalTN                        1.7861
## NodalStatusG1:LuminalLuminalB                1.4340
## NodalStatusG2:LuminalLuminalB                1.5342
## NodalStatusG3:LuminalLuminalB                0.7744
## NodalStatusG4:LuminalLuminalB                2.7401
## NodalStatusG1:LuminalTN                      1.8205
## NodalStatusG2:LuminalTN                      2.6375
## NodalStatusG3:LuminalTN                      1.5114
## NodalStatusG4:LuminalTN                      2.4145
## menostatPre:NodalStatusG1:LuminalLuminalB    9.5039
## menostatPre:NodalStatusG2:LuminalLuminalB    6.0287
## menostatPre:NodalStatusG3:LuminalLuminalB   13.3908
## menostatPre:NodalStatusG4:LuminalLuminalB        NA
## menostatPre:NodalStatusG1:LuminalTN          2.7586
## menostatPre:NodalStatusG2:LuminalTN          2.0828
## menostatPre:NodalStatusG3:LuminalTN          7.3393
## menostatPre:NodalStatusG4:LuminalTN        306.4633
## 
## Concordance= 0.596  (se = 0.012 )
## Rsquare= 0.067   (max possible= 1 )
## Likelihood ratio test= 54.69  on 28 df,   p=0.002
## Wald test            = 59.92  on 28 df,   p=4e-04
## Score (logrank) test = 69.53  on 28 df,   p=2e-05

Vemos que las combinaciones que mejor explcian el modelo son: menostat-tgrade, menostat-nodalStatus, tgrade-nodalStatus Haremos estas combinaciones teniendo como diferenciación principal el parámetro que ha demostrado ser el más significativo, menostat.

survfit(surv.obj ~ NodeStatus + menostat, BC, conf.type = "log-log") %>% 
    ggsurvplot(title = "Supervivencia estado nodal por menopausia", conf.int = T, 
        facet.by = "NodeStatus", legend.title = "Menostat", short.panel.labs = T)

survfit(surv.obj ~ NodalStatus + menostat, BC, conf.type = "log-log") %>% 
    ggsurvplot(title = "Supervivencia grupo nodal por menopausia", conf.int = T, 
        facet.by = "NodalStatus", legend.title = "Menostat", short.panel.labs = T)

survfit(surv.obj ~ Luminal+ menostat, BC, conf.type = "log-log") %>% 
    ggsurvplot(title = "Supervivencia Luminal por menopausia", conf.int=T, 
        facet.by = "Luminal", legend.title = "Menostat", short.panel.labs = T)

Y aplicando también la edad:

survfit(surv.obj~ NodeStatus+GrupoEdadCluster, BC, conf.type = "log-log") %>% 
    ggsurvplot(title = "Supervivencia de estado nodal por edad ", conf.int=T, palette = "Dark2", 
        facet.by = "GrupoEdadCluster")

survfit(surv.obj ~ NodalStatus+GrupoEdadCluster, BC, conf.type = "log-log") %>% 
    ggsurvplot(title = "Supervivencia entre Grupo nodal por edad", conf.int = F,  palette = c("slategray1", "#2E9FDF","blue2", "#E7B800","orange3"), facet.by = "GrupoEdadCluster", legend.title = "Grupo nodal", short.panel.labs = T)

survfit(surv.obj ~ Luminal +GrupoEdadCluster, BC, conf.type = "log-log") %>% 
    ggsurvplot(title = "Supervivencia entre Luminal por edad", conf.int=F, 
          facet.by = "GrupoEdadCluster", legend.title = "Luminal", short.panel.labs = T)

3.5.2.1 Supuestos observados

Diferenciado por menopausia:

  • Cuando la paciente es nodalmente positiva, ser premenopausica marca una diferencia en supervivencia ligeramente mayor, pero a simple vista es bastante similar que post.

  • Para el grado de los nodos: Como ya hemos visto tantas veces y seguimos ocmprobando, G3 marca la mayor diferencia para los niveles de menopausia, siendo la supervivencia de las Pre mucho mayor. También hay una diferencia fiable en G0. Para G4 encontramos destacable que es el único donde la supervivencia es mejor para Post.

  • Para el Luminal: Como hemos visto en su correspondiente tabla de riesgo, Para LuminalB ,para el tiempo 100), solo quedan tres pacientes vivas que viendo la tabla de riesgo sabemos que es en el grupo Post, por eso ese corte tan extraño en Pre. En cuanto al resto siguen una distribución dentro de lo esperado, para los 3 la esperanza de vida es claramente menor en Post.

En general para esta covariable la esperanza de supervivencia para las Pacientes Premenstruales es mejor en la mayoría de los casos, exceptuando Grupo nodal G4. Es especialmente mejor en Grupos nodales G3 y G0 y Nodalmente positivo.

Siguiendo esta misma linea, esperamos que las pacientes del grupo de edad más jóvenes tengas mejores expectativas, coincidiendo con las premenopausicas.

Diferenciando por edad:

  • Aquí encontramos la primera “contradicción”, así como veiamos que para las Pre la supervivencia era algo mejor con estado nodal positivo, viendo lor edad esto ocurre en el segundo grupo, donde todavia puede haber alguna menopausica, veamoslo:
with(BC, table(GrupoEdadCluster, menostat))
##                 menostat
## GrupoEdadCluster Post Pre
##          (25,50]   57 168
##          (50,66]  328  20
##          (66,90]  209   1

Efectivamente el número de pre en el segundo grupo de edad es mucho menor, pero no por ello inexistente, por lo que no deja de tener sentido el supuesto que habiamos hecho.

El grupo que cuenta con más pacientes ((50,66]) tiene más probabilidades de sobrevivir siendo Node negative, y sus intervalos de confianza son menores, por lo que en principio nos fiaremos más de este supuesto.

  • Grupo nodal podemos ver el patrón de como la curva decae antes para los grupos G3 y G4 que para el resto, coincidendo en menor o mayor medida para los 3 grupos de edad.

  • Luminal no hay tanta diferencia como en las anteriores variables. Podemos destacar que Triple Negativo tiene mejor supervivencia en el primer Grupo de edad.

Todos los resultados obtenidos parecen coherentes y cercanos a lo que esperábamos.


3.6. Pruebas de hipótesis del análisis exploratorio

Tendremos dos hipótesis:

En este caso, dado que hemos visto que los intervalos de confianza son amplios ya que algunos niveles cuentan con pocos datos, estableceremos un nivel de confianza del \(99%\).

  • \(H_0\) o hipótesis nula, la cual no rechazaremos a no ser que se demuestre lo contrario (\(p-value \leq 0,01\) o probabilidad del 99% de que \(H_0\) sea verdadera): Independencia de las curvas de supervivencia.

  • \(H_1\) o hipótesis alternativa, la cual aceptaremos en caso de rechazar \(H_0\): Dependencia de curvas.

La función survdiff() nos permite realizar comparación por parejas de forma no paramétrica de funciones de supervivencia, (Es la forma interna que tiene survminer de calcualar el \(pvalue\) que veiamos representado en las gráficas). Podremos observar sus correspondientes \(pvalue\) y \(\chi^2\).

Como hemos visto anteriormente, podemos usar para el ajuste de LogRank

  • rho=1: PetoPeto
  • rho=0: logRank

Consiste en realizar un test estadístico para ver si los resultados se deben al azar.

LogRank nos da la similitud entre dos grupos. La suma de las dos últimas columnas nos darán el valor de \(\chi^2\).

Chi Square: \(\chi^2=\frac{(0_1-E_1)^2}{E_1}+\frac{(0_2-E_2)^2}{E_2}\), explica como son os resultados observados en comparación con los esperados. Es otra manera de aceptar o rechazar la hipótesis nula

Estimación de riesgo: \(OR=\frac{01/E1}{02/E2}\)

survdiff(surv.obj~ NodeStatus + strata(menostat), BC, rho=0)
## Call:
## survdiff(formula = surv.obj ~ NodeStatus + strata(menostat), 
##     data = BC, rho = 0)
## 
##                            N Observed Expected (O-E)^2/E (O-E)^2/V
## NodeStatus=Node negative 347      325      356      2.72       5.6
## NodeStatus=Node positive 436      384      353      2.75       5.6
## 
##  Chisq= 5.6  on 1 degrees of freedom, p= 0.02
survdiff(surv.obj~ NodalStatus + strata(menostat), BC, rho=0)
## Call:
## survdiff(formula = surv.obj ~ NodalStatus + strata(menostat), 
##     data = BC, rho = 0)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## NodalStatus=G0 347      325    356.1     2.723     5.603
## NodalStatus=G1 108       96    102.8     0.451     0.534
## NodalStatus=G2 105       95     98.2     0.103     0.120
## NodalStatus=G3 185      163    129.5     8.678    10.894
## NodalStatus=G4  38       30     22.4     2.584     2.698
## 
##  Chisq= 15  on 4 degrees of freedom, p= 0.005
survdiff(surv.obj~ Luminal + strata(menostat), BC, rho=0)
## Call:
## survdiff(formula = surv.obj ~ Luminal + strata(menostat), data = BC, 
##     rho = 0)
## 
##                    N Observed Expected (O-E)^2/E (O-E)^2/V
## Luminal=LuminalA 538      493    501.9    0.1586    0.5480
## Luminal=LuminalB 110       98     86.2    1.6292    1.8800
## Luminal=TN       135      118    120.9    0.0707    0.0861
## 
##  Chisq= 1.9  on 2 degrees of freedom, p= 0.4
survdiff(surv.obj~ NodeStatus+strata(GrupoEdadCluster), BC,rho=0)
## Call:
## survdiff(formula = surv.obj ~ NodeStatus + strata(GrupoEdadCluster), 
##     data = BC, rho = 0)
## 
##                            N Observed Expected (O-E)^2/E (O-E)^2/V
## NodeStatus=Node negative 347      325      355      2.58      5.32
## NodeStatus=Node positive 436      384      354      2.59      5.32
## 
##  Chisq= 5.3  on 1 degrees of freedom, p= 0.02
survdiff(surv.obj~ NodalStatus+strata(GrupoEdadCluster), BC, rho=0)
## Call:
## survdiff(formula = surv.obj ~ NodalStatus + strata(GrupoEdadCluster), 
##     data = BC, rho = 0)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## NodalStatus=G0 347      325    355.3    2.5832    5.3188
## NodalStatus=G1 108       96    101.9    0.3455    0.4111
## NodalStatus=G2 105       95     97.5    0.0654    0.0768
## NodalStatus=G3 185      163    131.0    7.8125    9.8462
## NodalStatus=G4  38       30     23.2    1.9686    2.0617
## 
##  Chisq= 13.2  on 4 degrees of freedom, p= 0.01
survdiff(surv.obj~ Luminal +strata(GrupoEdadCluster), BC, rho=0)
## Call:
## survdiff(formula = surv.obj ~ Luminal + strata(GrupoEdadCluster), 
##     data = BC, rho = 0)
## 
##                    N Observed Expected (O-E)^2/E (O-E)^2/V
## Luminal=LuminalA 538      493      504   0.22650    0.7970
## Luminal=LuminalB 110       98       88   1.14202    1.3298
## Luminal=TN       135      118      117   0.00368    0.0045
## 
##  Chisq= 1.4  on 2 degrees of freedom, p= 0.5
Supuesto Chisq pvalue Hipótesis
NodalStatus + menostat 15 0.0005 \(H_1\)
NodalStatus+GrupoEdadCluster 13.2 0.01 \(H_1\)
NodeStatus + menostat 5.6 0.02 \(H_0\)
NodeStatus+GrupoEdadCluster 5.3 0.02 \(H_0\)
Luminal + menostat 1.9 0.4 \(H_0\)
Luminal +GrupoEdadCluster 1.4 0.5 \(H_0\)

Aceptamos \(H_0\) (independencia de las curvas) para un 99% de confianza para todas excepto para las estratificacionesmenopausia y edad para estado nodal, las cuales han rechazalo la hipótesis nula por tener un \(p-value<0,01\) y aceptamos su dependencia, entendiendo que no dan tanta información sobre el riesgo.

Parece que la estratificación por grupo de edad marca más diferencia que la menopausia, son grupos más específicos.

Nos quedamos con que la mejor explicación del riesgo viene dada por Luminal y Estado nodal (Positivo/Negativo), todos estratificados por los grupos de edad obtenidos de la clusterización.

survfit(surv.obj ~ Luminal +GrupoEdadCluster, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia entre grado tumoral por edad", conf.int=F, facet.by = "GrupoEdadCluster", legend.title = "Grupo nodal", short.panel.labs = T, fun="cumhaz")

survfit(surv.obj~ NodeStatus+GrupoEdadCluster, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia de estado nodal por edad ", conf.int=F, palette = "Dark2", facet.by = "GrupoEdadCluster", fun="cumhaz")


3.7. Análisis de estratificación de edad

3.7.1. Parametricamente vs No parametricamente

levels(BC$GrupoEdadCluster)
## [1] "(25,50]" "(50,66]" "(66,90]"
grupoedad1<-BC %>% filter(GrupoEdadCluster=="(25,50]")
grupoedad2<-BC %>% filter(GrupoEdadCluster=="(50,66]")
grupoedad3<-BC %>% filter(GrupoEdadCluster=="(50,66]")
grupos=list(grupoedad1,grupoedad2,grupoedad3)

params<- lapply(grupos, function(x) flexsurvreg(surv.obj ~1, data=x, dist = "exp") %>% summary(type="survival") %>% data.frame)

no.params<-lapply(grupos, function(x) survfit(surv.obj~ 1, data=x) %>% fortify)

gglist<-list()
gglist.risk<-list()
titles=c("(25,50]","(50,66]","(50,66]")
for(model in 1:length(grupos)) {
  gglist[[model]]<-ggplot() + geom_line(aes(time, est, col = "Parametrico"), data=params[[model]]) + geom_step(aes(time, surv, col = "No Parametrico"), data=no.params[[model]]) + labs(x = "Tiempo", y = "Probabilidad de Supervivencia", col = "Ajustes", title = titles[model])
  
  gglist.risk[[model]]<-ggplot() + geom_line(aes(time, -log(est), col = "Parametrico"), data=params[[model]]) + geom_step(aes(time, -log(surv), col = "No Parametrico"), data=no.params[[model]]) + labs(x = "Tiempo", y = "Riesgo acumulado", col = "Ajustes", title = titles[model])
}
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
grid.arrange(gglist[[1]], gglist[[2]],gglist[[3]], ncol=2, nrow=2)

grid.arrange(gglist.risk[[1]], gglist.risk[[2]],gglist.risk[[3]], ncol=2, nrow=2)

3.7.2. Estadísticamente

Los parámetros AIC, BIC y LogLik ya están explicados en el punto \(3.3.\)

params<- sapply(grupos, function(x) flexsurvreg(surv.obj ~1, data=x, dist = "exp"), USE.NAMES = T, simplify = F)
no.params<- sapply(grupos, function(x) survfit(surv.obj ~1, data=x), USE.NAMES = T, simplify = F)

means.param <- sapply(params, function(x) print(x, print.rmean=T))
## Call:
## flexsurvreg(formula = surv.obj ~ 1, data = x, dist = "exp")
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  0.024097  0.022387  0.025937  0.000905
## 
## N = 783,  Events: 709,  Censored: 74
## Total time at risk: 29422.92
## Log-likelihood = -3350.503, df = 1
## AIC = 6703.005
## 
## Call:
## flexsurvreg(formula = surv.obj ~ 1, data = x, dist = "exp")
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  0.024097  0.022387  0.025937  0.000905
## 
## N = 783,  Events: 709,  Censored: 74
## Total time at risk: 29422.92
## Log-likelihood = -3350.503, df = 1
## AIC = 6703.005
## 
## Call:
## flexsurvreg(formula = surv.obj ~ 1, data = x, dist = "exp")
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  0.024097  0.022387  0.025937  0.000905
## 
## N = 783,  Events: 709,  Censored: 74
## Total time at risk: 29422.92
## Log-likelihood = -3350.503, df = 1
## AIC = 6703.005
means.noparam<- sapply(no.params, function(x) print(x, print.rmean=T))
## Call: survfit(formula = surv.obj ~ 1, data = x)
## 
##          n     events     *rmean *se(rmean)     median    0.95LCL 
##      783.0      709.0       40.7        1.3       27.6       24.5 
##    0.95UCL 
##       31.0 
##     * restricted mean with upper limit =  217 
## Call: survfit(formula = surv.obj ~ 1, data = x)
## 
##          n     events     *rmean *se(rmean)     median    0.95LCL 
##      783.0      709.0       40.7        1.3       27.6       24.5 
##    0.95UCL 
##       31.0 
##     * restricted mean with upper limit =  217 
## Call: survfit(formula = surv.obj ~ 1, data = x)
## 
##          n     events     *rmean *se(rmean)     median    0.95LCL 
##      783.0      709.0       40.7        1.3       27.6       24.5 
##    0.95UCL 
##       31.0 
##     * restricted mean with upper limit =  217

La diferencia no es apreciable en ningún parámetro.

3.7.3. Cox

Proporcionalidad de los riesgos

cox.edad<-coxph(surv.obj ~ GrupoEdadCluster, BC)
check.ph<-cox.zph(cox.edad)
check.ph
##                             rho chisq      p
## GrupoEdadCluster(50,66] -0.0189 0.251 0.6162
## GrupoEdadCluster(66,90] -0.0747 3.862 0.0494
## GLOBAL                       NA 4.276 0.1179
ggcoxzph(cox.zph(cox.edad))

Volvemos a comprobar que la hipótesis de riesgos proporcionales es válida para la estratificación de edad, y visualizamos la distrbución.

Schoenfeld residuals:

plot(check.ph, var = 1)
abline(h = coef(cox.edad)[1], col = "red", lwd = 2)

plot(check.ph, var = 1)
abline(h = coef(cox.edad)[2], col = "red", lwd = 2)

Vemos los 3 grupos claramente diferenciados, solo un de ellos encajando.

Análisis Residuales

ggcoxdiagnostics(cox.edad, type = "dfbeta")

Podemos ver que bastantes punyos tienen mucha influencia en la linea de regresión, y también tenemos algunos outlayers (más para el grupo de edad mayor)

ggcoxdiagnostics(res.cox, type = "deviance", ox.scale = "time")

4.Conclusiones